REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202^4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503. 


1.  AQmCy  \JSEOHV<  (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

25Jun.02  DISSERTATION  _ 


4.  TITLE  AND  SUBTITLE  5.  FUNDING  NUMBERS 

SHEAR  LAYER  STABILITY  IN  A  TWO  DIMENSIONAL  DISK 


6.  AUTHOR(S) 

CAPT  WOLVERTON  ROBERT  H 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

UNIVERSITY  OF  NEW  MEXICO 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


CI02-147 


9.  SPONSORING/MONITORING  AGENCY  NAME{S)  AND  ADDRESS(ES) 

THE  DEPARTMENT  OF  THE  AIR  FORCE 

AFIT/CIA,  BLDG  125 

2950  P  STREET 

WPAFB  OH  45433 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


12a.  DISTRIBUTION  AVAILABILITY  STATEMENT 

Unlimited  distribution 

In  Accordance  With  AFI  35-205/AFIT  Sup  1 


12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  (Maximum  200  words) 


20020702  031 


15.  NUMBER  OF  PAGES 

113 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRAC 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 


Standard  Form  298  (Rev.  2-89  (EG) 

Prescribed  by  ANSI  Std.  239.18 

Designed  using  Perform  Pro,  WHS/DIOR,  Oct  94 


SHEAR  LAYER  STABILITY 
IN  A  TWO-DIMENSIONAL  DISK 


BY 

ROBERT  H.  WOLVERTON 

B.S.,  Engineering  Sciences,  USAF  Academy,  1991 
M.S.,  Aeronautical  Engineering,  University  of  Washington,  1992 


DISSERTATION 

Submitted  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 

Doctor  of  Philosophy  in  Mathematics 


The  University  of  New  Mexico 
Albuquerque,  New  Mexico 


December,  2001 


SHEAR  LAYER  STABILITY 
IN  A  TWO-DIMENSIONAL  DISK 


ROBERT  H.  WOLVERTON 


Doctor  of  Philosophy  in  Mathematics 


December,  2001 


11 


DEDICATION 


In  memory  of  my  brother,  Paul  F.  Wolverton. 
His  abilities  far  exceeded  my  own. 


iii 


ACKNOWLEDGEMENTS 


My  heartfelt  and  deepest  appreciation  of  Professor  Evangelos  A.  Coutsias  who 
provided  everything  required  for  this  journey  from  faith  and  guidance  to  grit  and 
determination.  He  is  without  peer. 

I  am  also  grateful  for  my  committee  members  and  their  service.  Thanks  to  Prof. 
Mclver,  Prof.  Hagstrom,  Prof.  Kapitula,  Dr.  Romero,  and  Dr.  Alsing.  Furthermore,  I’d 
like  to  thank  Dr.  Nielsen  for  his  support. 

All  the  best  plans  require  implementation.  To  that  end  my  sincerest  thanks  to  Ms. 
Roxanne  Littlefield;  she  can  bring  any  bureaucracy  to  its  knees  with  a  single  glance.  And 
a  brass  ovation  for  the  only  men  who  could  solve  my  perpetual  computer  dilemma,  Mr. 
Dann  “Thaiday”  Brewer  and  Mr.  Craig  “Duct  Tape”  Lewis,  the  most  resourceful  sys¬ 
admins  in  town. 

Thanks  to  USAFA  Dept  of  Mathematical  Sciences  for  providing  me  the 
opportunity  to  pursue  my  dream.  This  work  was  partially  supported  by  the  National 
Computational  Science  Alliance  under  DMS000002N  using  the  alliance  AHPCC  Linux 
cluster  and  by  the  National  Science  Foundation  under  DMS9977396. 


IV 


SHEAR  LAYER  STABILITY 
IN  A  TWO-DIMENSINAL  DISK 


BY 

ROBERT  H.  WOLVERTON 


ABSTRACT  OF  DISSERTATION 

Submitted  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 

Doctor  of  Philosophy  in  Mathematics 


The  University  of  New  Mexico 
Albuquerque,  New  Mexico 


December,  2001 


TABLE  OF  CONTENTS 


LIST  OF  FIGURES . ix 

1  Introduction .  1 

1.1  Kelvin-Helmholtz  instability .  5 

1.2  Problem  formulation .  9 

1.2.1  Equations  of  motion .  10 

1.2.2  Ekman  friction .  15 

1.2.3  Energy  evolution .  16 

1.2.4  Enstrophy  evolution .  19 

1.2.5  Angular  momentum .  21 

1.2.6  Circulation  evolution .  23 

1.2.7  Circulation  evolution .  25 

2  Numerics . 27 

2.1  Introduction .  27 

2.2  Fourier-Chebyshev  expansion . 29 

2.3  Spatial  differentiation  in  the  spectral  scheme .  31 

2.4  Postconditioning .  33 

2.5  Computational  algorithm .  35 

2.5.1  Initial  conditions .  36 

2.5.2  Forcing  of  the  flow .  37 

2.5.3  Boundary  conditions .  39 

2.5.4  Accurate  initiation  of  the  time  scheme . 41 

2.5.5  Time  integration  scheme . 42 

2.5.6  Diagnostics . 43 

2.6  Forced  periodicity . 44 

2.7  Conditioning  of  the  problem . 44 

3  Stability  analysis . 46 

3.1  Introduction . 46 

vii 


47 


3.2  Homogeneous  flow  . 

3.3  Solving  the  generalized  eigenvalue  problem,  .  52 

3.4  Secondary  bifurcations .  54 

3.5  Solving  the  generalized  eigenvalue  problem.  Ad),  =  .  59 

4  Results . 61 

4.1  Introduction . 61 

4.2  Homogeneous  instability . 62 

4.3  Asymmetric  tripole .  69 

4.4  Secondary  instability . 74 

4.4. 1  Period  doubling .  74 

4.4.2  Oscillatory  states .  82 

APPENDIX  A:  Vector  identities .  92 

APPENDIX  B:  Conditioning  of  the  problem .  93 

APPENDIX  C:  Chebyshev  recurrence  relations .  97 

APPENDIX  D:  Cheby she V  matrices .  99 

APPENDIX  E:  Chebyshev  convolution  matrix . 102 

APPENDIX  F:  Fourier  convolution  matrix . 104 

BIBLIOGRAPHY . Ill 


viii 


LISTOFHGURES 


1.1  Disk  geometry .  4 

1.2  Kelvin-Helmholtz  instability .  5 

1.3  Stewartson  layers . 12 

2. 1  Velocity  at  the  shear  layer . 38 

2.2  Vorticity  at  the  shear  layer . 38 

3.1  Perturbation  analysis . 47 

3.2  Homogeneous  pseudospectrum . 53 

4. 1  Shear  layer  in  the  annulus . 62 

4.2  Disk  v.s.  annulus  marginal  stability,  F  =  6.5 . 64 

4.3  Disk  v.s.  annulus  marginal  stability,  F  =  10 . 64 

4.4  Marginal  stability  . 65 

4.5  Spectrum  for  Hopf  bifurcation . 66 

4.6  Marginal  stability  with  Ekman  comparison . 67 

4.7  Flow  after  Hopf  bifurcation.  Re  =  40 . 68 

4.8  Flow  of  a  point  vortex.  Re  =  250  . 70 

4.9  Energy  and  enstrophy  for  a  point  vortex.  Re  =  250  . 71 

4. 10  Flow  of  a  point  vortex.  Re  =  150 . 72 

4. 1 1  Energy  and  enstrophy  for  the  transition  to  a  point  vortex . 73 

4.12  Spectrum  of  a  inhomogeneous  flow.  Re  =  40 . 76 

4.13  Period  doubling  bifurcation,  n  =  6  —>  n  =  3 . 77 

4.14  Energy  and  enstrophy  for  the  n  =  6  — >  «  =  3  bifurcation . 78 

4.15  Re  =  75,  n  =  5  braid . 79 

4.16  Energy  and  enstrophy  for  initial  transition  to  n  =  5  braid . 80 

4.17  Re  =  75 ,  forced  periodicity  n  =  6  braid . 81 

4.18  Bifurcation  steps  from  n  =  4  — >n  =  3 . 83 

4.19  4  ®  2  periodic  flow . 84 

4.20  Energy  and  enstrophy  for  the  4  ©  2  periodic  flow . 85 

4.21  4  ©  3  ©  2  periodic  flow . 86 

4.22  Energy  and  enstrophy  for  the  4  ©  3  ©  2  periodic  flow . 87 


IX 


4.23  Still  frame  showing  the  mode  3  impact  on  4  ©  3  ®  2  periodic  flow . 88 

4.24  Energy  and  enstrophy  for  the  transition  from  4  ©  2  steady  to 

4® 3® 2  periodic . 89 

4.25  Irregular  flow . 90 

4.26  Energy  and  enstrophy  for  the  irregular  flow . 91 


X 


SHEAR  LAYER  STABILITY 
IN  A  TWO-DIMENSIONAL  DISK 


by 

R.  H.  Wolverton 

B.S.,  Engineering  Sciences,  US  Air  Force  Academy,  1991 
M.S.,  Aeronautical  Engineering,  University  of  Washington,  1992 
Ph.D.,  Mathematics,  University  of  New  Mexico,  2001 

ABSTRACT 

The  dynamics  of  rotating  fluids  provide  a  rich  aggregate  of  periodic,  quasi- 
periodic,  and  irregular  behavior.  Many  investigations  of  two-dimensional  (2D)  flows 
containing  fluid  velocity  inflections  present  the  Kelvin-Helmholtz  (KH)  instability.  In 
this  investigation  we  study  the  saturation  of  the  KH  instability  for  a  forced  circular  shear 
layer  in  a  differentially  rotating  split-disk. 

Complex  vortex  interactions  are  reasonably  well  understood  through 
experimentation  but  modeling  them  requires  highly  accurate  numerical  schemes.  To 
explore  these  flows  our  investigations  employ  an  efficiently  parallelized,  highly  accurate 
pseudospectral  scheme  for  the  solution  of  the  incompressible  Navier-Stokes  equations  in 
a  disk  geometry,  Torres  &  Coutsias  [32].  Beyond  the  initial  KH  instability,  secondary 
transitions  in  the  flow  yield  symmetry  breaking  bifurcations  resulting  in  periodic  and 
irregular  states.  Simulations  are  provided  for  the  intermediate  states  between  the  n  =  4 
and  n  =  3  vortex  braids.  Oscillating  states  not  previously  seen  in  numerical  studies  are 
reported.  Unlike  the  jump  transitions  between  braids  of  different  order,  the  oscillating 
states  were  found  to  be  supercritical  bifurcations  and  thus,  not  hysteretic.  Period 
doubling  bifurcations  are  observed  during  some  spin-up  studies  in  which  intermediate 
symmetry  breaking  bifurcations  are  bypassed. 

Accurate  spectral  simulation  offers  the  means  for  systematic  exploration  of  the 
dynamics  associated  with  rotating  fluids.  Herein  we  construct  such  a  scheme  and  present 
bifurcation  analysis  for  secondary  transitions. 
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1  Introduction 


Rotating  shear  layers  occur  in  a  variety  of  situations  from  astronomical 
phenomena  to  industrial  applications.  The  Great  Red  Spot  on  Jupiter  and  the  Blue  Spot 
on  Neptune  result  from  differential  rotation  associated  with  sharply  sheared  zonal  flows 
within  their  respective  atmospheres.  Similar  vortices  are  created  in  Earth’s  atmosphere 
through  the  instabilities  of  differential  velocity  associated  with  fluid  propagation.  Schar 
&  Davies  [28]  and  Simmons  &  Hoskins  [29]  considered  weather  system  evolution  and 
baroclinic  front  instabilities.  Weber  &  Smith  [34]  followed  the  life  cycle  of  tropical 
cyclones,  and  additional  work  by  Read  et  al  [27]  considers  these  instabilities  combined 
with  a  temperature  gradient.  Turbo-machinery  within  a  closed  housing  produces  stresses 
via  instabilities  in  the  shear  layer  created  by  the  high  speed  rotation.  One  example  is 
computer  hard  drives,  and  Humphrey  &  Gore  [19]  consider  the  design  of  these  drives  and 
the  effects  of  shear  layer  instabilities. 

In  planar  and  circular  2-D  flow  the  Kelvin-Helmholtz  instability  arises  at  a  shear 
layer  with  an  inflection  in  the  fluid’s  velocity  profile.  A  derivation  of  the  instability  is 
provided  in  section  1.1.  The  dynamics  associated  with  the  Kelvin-Helmholtz  instability 
apply  to  both  viscid  and  inviscid  fluids  with  viscosity  adding  disspative  capacity  to  the 
fluid  as  interned  friction.  Nonlinear  evolution  of  a  circular  shear  layer  in  a  rotating  fluid 
has  been  investigated  in  several  forms  beginning  with  Stewartson’s  [30]  analysis  of  a 
small  differential  rotation  in  a  strong  background  flow.  Greenspan  [17]  expanded 
Stewartson’s  work  by  presenting  in-depth  analysis  of  rotating  fluids.  Of  particular 
interest  are  the  zones  of  vertical  fluid  transport,  the  Stewartson  layers,  that  develop  near 
the  shear  layer  in  a  differentially  rotating  viscous  flow. 

Experimental  verification  of  shear  layer  instabilities  has  been  demonstrated  in 
several  geometries.  While  Tomasini  et  al  [31]  consider  the  effects  of  compressibility,  all 
the  work  discussed  throughout  this  paper  is  conducted  using  low  speed,  incompressible 
fluids.  Hide  &  Titman  [18]  completed  experiments  in  a  rotating  tank  equipped  with  a 
differentially  rotating  disk.  Their  experiment  showed  a  well  defined  parameter  value  for 
the  onset  of  stability,  and  a  vorticity  mode  number  that  decreases  as  the  differential 
velocity  in  the  disk  is  increased.  Flow  within  a  thin  annual  region  was  investigated  by 
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Rabaud  &  Couder  [25]  who  formulated  a  correlation  between  the  initial  instability  and 
the  Reynolds  number  of  the  flow.  Their  work  was  extended  by  Chomaz  et  al  [6]  giving  a 
deeper  analysis  of  the  secondary  transitions  occurring  after  the  initial  bifurcation.  With  an 
experimental  apparatus  similar  to  the  one  used  by  Hide  &  Titman  [18],  Niino  &  Misawa 
[23]  characterized  the  initial  nonlinear  stability  of  the  flow  including  Stewartson  layers. 
The  effects  of  both  internal  viscous  dissipation  and  Ekman  pumping  are  analyzed 
providing  insights  about  the  flow  within  the  Stewartson  layer.  Konijnenberg  et  al  [20] 
continue  the  experimental  investigations  using  a  parabolic  bowl.  This  non-annular,  free 
surface  parabolic  bowl  allows  for  the  formation  of  a  point  vortex  flow.  In  a  near  disk 
geometry,  Friih  &  Read  [15]  experimentally  analyze  secondary  transitions  using  detailed 
flow  visualization  of  large  scale  structures  via  laser-Doppler  velocimetry.  Experiments 
completed  by  Hide  &  Titman  [18],  Niino  &  Misawa  [23],  Konijnenberg  et  al  [20],  and 
Friih  &  Read  [15]  use  rotating  tanks  where  the  fluid  layer  is  thicker  than  the  Ekman  layer 
thickness.  In  contrast  Rabaud  &  Couder  [25]  and  Chomaz  et  al  [6]  impose  a  restriction 
on  the  fluid  thickness  as  less  than  the  Ekman  layer  thickness. 

Recent  studies  have  investigated  quasi-two-dimensional  (Q2D)  flow  which 
includes  the  effects  of  bottom  friction  and  the  resulting  Ekman  layer.  Problems  are 
classified  as  Q2D  when  the  vertical  fluid  velocity  is  represented  by  an  approximation  of 
the  planar  velocity.  Danilov  et  al  [9]  validate  the  Q2D  approximation  for  Ekman  and 
Rayleigh  friction  by  comparing  experimental  and  analytical  results.  The  Stewartson 
layers  which  support  vertical  fluid  transport  near  the  shear  layer  have  been  investigated 
by  several  authors.  Hide  &  Titman  [18],  Niino  &  Misawa  [23],  Konijnenberg  et  al  [20], 
etc.  Section  1.2.2  provides  a  detailed  discussion  of  the  Stewartson  layers.  Energy 
cascades  studied  by  Manin  [21]  based  on  Komolgorov  wave  numbers  give  a  correlation 
between  the  wave  number  and  energy  cascade  rates.  Inverse  energy  cascade  and  self¬ 
organization  dominate  the  dynamics  of  a  rotating  fluid  in  a  bounded  domain.  Churilov  & 
Shukhman  [7]  analyze  and  classify  the  initial  bifurcation  from  steady,  homogeneous 
rotational  flow.  Bergeron  et  al  [3]  use  weakly  nonlinear  stability  analysis  to  validate 
numerical  simulation  of  the  initial  bifurcation.  A  perturbation  of  the  initial  flow  is 
applied  to  the  azimuthally  symmetric  state.  The  resulting  perturbed  problem  is  solved  in 
its  linear  form  giving  the  critical  Reynolds  number  for  the  bifurcation  from 
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homogeneous,  steady  flow.  The  third  order  Landau  equations  of  the  weakly  nonlinear 
theory  yield  the  form  of  the  bifurcation  as  a  supercritical  pitchfork.  A  linear  analysis 
ignoring  viscosity  conducted  by  Konijnenberg  et  al  [20]  considers  the  effect  on  the 
initial  instability  of  a  rotating  fluid. 

Secondary  bifurcations  from  a  steady  n  vortex  flow  have  been  the  focus  of 
several  experiments  and  numerical  simulations.  Chomaz  et  al  [6]  following  the  work  of 
Rabaud  &  Couder  [25]  determine  two  distinctly  different  forms  of  instability  depending 
on  the  initial  number  of  vortices.  They  explain  two  examples,  n-5  and  n  =  4 ,  in  detail. 
Transition  from  a  solution  composed  of  an  even  number  of  vortices  begins  with  the  first 
subharmonic  symmetry  breaking.  For  a  solution  composed  of  an  odd  number  of  vortices 
the  first  subharmonic  is  not  discretized  along  the  shear  layer  so  a  traveling  wave 
develops.  In  general  the  amplification  of  subharmonic  and  superharmonic  waves  precede 
the  secondary  transitions.  The  numerical  analysis  by  Rabaud  &  Couder  [25]  and  Chomaz 
et  al  [6]  uses  a  double  Fourier  grid  at  low  resolution  and  only  captures  large  scale 
transitions.  Supporting  the  Chomaz  et  al  [6]  work,  Bergeron  et  al  [3]  use  a  Fourier- 
Chebyshev  expansion  to  provide  accurate,  spectral  numerical  simulation  of  the  initial 
bifurcation  in  an  annular  geometry.  Investigation  by  Konijnenberg  et  al  [20],  also  using  a 
Fourier-Chebyshev  expansion  in  an  annular  domain,  yields  simulations  of  the  p  effect 
that  match  their  experimental  findings.  Niino  &  Misawa  [23]  study  the  barotropic 
instability  and  find  supercriticality  of  the  flow  and  a  decrease  in  the  wavenumber  of  the 
solution.  However,  their  analysis  predicts  increasing  and  decreasing  wavenumbers  as  the 
Reynolds  number  of  the  flow  is  increased.  Danilov  et  al  [9]  use  a  double  Fourier 
expansion  to  simulate  magnetohydrodynamic  flow  in  a  rectangular  domain  with  emphasis 
given  to  vortex  interactions.  Friih  &  Read  [15]  parameterize  their  experiment  for  a  wide 
range  of  Rossby  and  Ekman  numbers.  Using  the  parameterization  particular  attention  is 
given  to  the  secondary  bifurcation  regions  and  period  doubling.  The  vortices  are  mapped 
via  particle  tracking  giving  a  detailed  map  for  the  evolution  of  a  vortex.  Manin  [21] 
analyzes  the  merger  of  two  vortices  providing  insight  into  the  transition  to  lower 
frequency  solutions  within  rotating  fluids.  Rand  [26]  described  the  secondary  transitions 
via  symmetry  breaking  bifurcations.  He  reports  that  during  the  transition  the  group 
dimension  of  initially  continuous  spatio-temporal  symmetries  drops  by  one  in  essence 
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collapsing  to  a  non-trivial  lattice  in  the  plane.  Flow  evolution  through  the  transitions  is 
qualified  in  the  following  progression:  homogeneous,  periodic,  quasi-periodic,  and 
turbulent. 

The  Kelvin-Helmholtz  instability  of  a  viscous,  incompressible  fluid  is  a  the 
starting  point  for  a  series  of  rich  bifurcations  in  the  2D  circular  geometry.  Our  circular 
geometry  is  not  annular  but  rather  is  a  split  disk  with  forcing  introduced  through  Ekman 
pumping,  see  figure  1.1. 


h 


fig  1.1  The  disk  geometry  forced  with  Ekman  pumping  from  the  top  and  bottom. 
The  shear  layer  is  forced  at  a  distance  a  and  the  height  of  the  cell  is  h. 


The  Ekman  pumping  is  induced  by  rotating  the  inner  portion  of  the  top  and 
bottom  walls.  The  flow  is  characterized  as  Q2D  where  the  vertical  fluid  transport  due  to 
the  Ekman  forcing  is  included  with  a  planar  approximation.  Stability  of  the  initial 
bifurcation  is  resolved  demonstrating  the  similarity  between  the  disk  and  annulus 
geometries  for  low  Reynolds  numbers,  and  verification  of  the  results  is  demonstrated 
both  analytically,  Bergeron  et  al  [3],  and  experimentally,  Fruh  &  Read  [15].  Secondary 
transitions  are  investigated  with  emphasis  on  symmetry  breaking  bifurcations  (including 
period  doubling  bifurcations)  and  periodic  states.  Using  an  accurate,  pseudospectral 
scheme,  we  provide  the  first  simulation  of  the  4  ©  3  ®  2  periodic  state  and  its  associated 
transitions.  The  notation  4  ©  3  ©  2  implies  the  flow  is  combination  of  Fourier  modes  4,  3, 
and  2.  The  4  ©  3  ©  2  state  consists  of  temporally  periodic  flow  built  on  the  interaction  of 
the  n  =  4  ©  2  and  n  =  3  solutions.  Also,  we  provide  the  first  numerical  investigation  of 
irregular  flow  regimes.  The  irregular  flow  persists  for  a  narrow  range  of  Reynolds 
number  beyond  the  4  ©  3  ©  2  periodic  state.  Finally,  we  simulate  the  asymmetric  tripole, 
and  introduce  the  complex  nature  of  its  parameter  space.  The  sections  of  this  thesis 
include  (2)  Numerics,  (3)  Stability  analysis,  and  (4)  Results. 
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1.1  Kelvin-Helmholtz  instability 

This  phenomenon  results  from  the  intrinsic  instability  of  a  sheet  vortex  to  small 
perturbations.  Transition  layer  instability  results  in  vortex  roll-up  to  the  well  documented 
vortex  pairing,  e.g.  Winant  &  Browand  [35].  The  Kelvin-Helmholtz  instability  is 
characterized  by  an  inflection  of  the  fluid  velocity  at  the  critical  layer,  Chandrasekhar  [5]. 
The  inflection  creates  increasing  undulations  at  the  shear  layer  driving  the  flow  to 
turbulence  in  the  3D  scenario.  An  analysis  of  the  instability  is  presented  here  to  provide 
background  information  on  a  dynamical  instability  possessed  by  many  unidirectional 
flows. 

A  simple  model  builds  upon  Bachelor  [2]  to  demonstrate  the  instability  has 
irrotational  flow  on  both  sides  of  a  slightly  perturbed  vortex  sheet.  For  consistency  with 
the  work  presented  in  this  paper,  the  flow  is  considered  within  a  polar  coordinate  system. 
The  Kelvin-Helmholtz  instability  is  a  local  phenomenon,  and  it  suffices  to  consider  a  thin 
sector  of  the  2-dimensional  disk  and  focus  on  the  area  near  the  shear  layer.  While  the 
shear  layer  in  the  disk  geometry  has  a  nontrivial  curvature,  the  shear  layer  in  a  thin  slice 
of  the  flow  is  nearly  linear.  The  flow  is  considered  in  two  regions  around  the  vortex 
sheet,  outer  and  inner  see  figure  1.2. 


r 


fig  1.2  A  perturbed  vortex  sheet  near  the  shear  layer.  7]  is 
the  distance  between  the  vortex  sheet  and  the  ideal  shear  layer. 


Since  the  flow  is  irrotational  and  incompressible,  the  velocity  is  written  in  terms 
of  a  potential  function.  The  outer  flow  has  a  velocity,  +  V^j ,  and  the 

inner  flow  has  a  velocity,  +  V^2 »  where  the  disturbance  potentials,  tp^  and 

(P2 ,  are  small.  The  distance  between  the  slightly  deformed  vortex  sheet  and  its  ideal 
position  constitutes  a  small  perturbation  varying  with  the  angle,  6 ,  and  with  time 
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All  of  the  small  quantities,  rj,  (p^,  and  <P2  are  related  through  the  material  boundary 
created  by  the  vortex  sheet.  A  few  items  to  note  before  continuing  are: 

1)  the  gradient  of  the  perturbation  is 

2)  the  total  derivative  with  respect  to  time  is 

3)  the  boundary  formed  by  the  vortex  sheet  is  located  at  a +  7j  where  a  is  the 
distance  from  the  center  of  the  disk  to  the  ideal  shear  layer.  Any  fluid  parameter 
^  evaluated  at  the  boundary  is  written  in  the  form 

and  for  the  first  order  approximations  the  magnitude  of  the  perturbation  is  nearly 
zero  and  the  notation  for  the  boundary  is  defined  as 

We  begin  by  noting  the  relationship  between  the  perturbation  velocity  and  the 
gradient  of  the  disturbance  potentials  in  the  radial  direction.  Velocity  in  the  radial 
directions  results  from  fluctuations  in  the  perturbation  and  for  the  outer  flow  this  gives 

L  =  ^  =  )|^  ■  V;/ 

which  simplifies  after  completing  the  dot  product 

Considering  the  disturbances  and  the  perturbation  to  be  small  parameters,  the  first  order 
approximation  for  the  outer  flow  velocity  gradient  in  the  vertical  direction  yields 

(1.1)  d,(Pil^=d,t7-^U0d0?j 

Similarly  the  radial  velocity  of  the  inner  flow  resulting  from  the  evolution  of  the  vortex 
sheet  perturbation  is 

(1.2)  d,<Pil^=d,rj  +  -^Ugd0f? 

In  addition  to  the  material  connectivity  at  the  vortex  sheet,  the  pressure  gradient  across 
the  shear  layer  is  zero.  As  a  result  the  difference  in  pressure  on  either  side  of  the  vortex 
sheet  is  zero.  For  the  incompressible,  irrotational  flow  pressure  is  given  according  to 
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(1.3)  p  =  pQ-p{d,(p  +  \u-u) 

Combining  this  definition  for  pressure  with  the  homogeneity  of  the  flow  yields  continuity 
for  the  velocity  potential  evolution  rate  across  the  shear  layer.  Here  is  the  result  via  the 
pressure  difference 

(Pl  )|g  ®  " ^ou(«r Ig 

and  expanding  the  velocity  products  gives 

+\uede(p2  +j^{^6(P2f\  =0 

The  first  order  approximation  for  the  above  equation  leaves  only  the  relationship  between 
the  time  partials  for  the  disturbance  potentials  and  their  associated  azimuthal  gradients. 
No  radial  gradients  are  present  in  the  first  order  rendering 

-4)  la  “  ^t^Pi  la  “  2a  15  “  '^^0^2 1 j  “  ® 

The  equations  (1.1),  (1.2),  and  (1.4)  can  be  solved  if  the  perturbation  and  disturbances  are 
expanded  with  a  Fourier  basis.  It  is  evident  that  a  sinusoidal  dependence  of  ;/  on  0 
demands  a  similar  dependence  from  the  disturbance  potentials.  Thus  all  the  small 
parameters  are  represented  by  the  superposition  of  Fourier  modes 

77,^1, ^2  “= 

where  each  k  mode  behaves  independently.  The  disturbance  potentials  are  restricted  to  a 
subset  of  the  solutions  for  (1.1),  (1.2),  and  (1.4)  through  the  irrotational  nature  of  the 
fluid.  Namely  (p^  and  ^2  must  independently  satisfy  Laplace’s  equation 

V  V;  =  r^dj(Pj  +  rd,(Pj  +  dl<Pj  =  0 

where  7  =  1, 2 .  In  general  Laplace’s  equation  in  polar  coordinates  can  be  solved  using  a 
product  solution  while  considering  the  radial  and  azimuthal  components  separately.  The 
solutions  are 

«;,  =  (c/+c,r-‘y“ 
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The  solutions  for  (p^  and  ^2  are  consistent  with  the  Fourier  modal  assumption.  To 
simplify  these  solutions,  we  enforce  conditions  such  that  the  disturbance  motion  vanishes 
away  from  the  vortex  sheet.  Thus  the  outer  flow  disturbance  potential  vanishes  for 
increasing  values  of  r  implying  Cj  =  0 .  For  the  inner  flow  the  radius  approaches  zero 
requiring  d2=0.  After  establishing  the  disturbance  potentials  as  solutions  for  Laplace’s 
equation,  the  time  dependence  of  the  potentials  is  reincorporated.  Writing  the  solutions 
using  time  dependent  coefficients  produces 

(1.5)  * 

The  perturbation,  r] ,  is  not  a  function  of  the  radius  and  when  expanded  in  Fourier  modes 
takes  on  the  form 

(1.6)  7]  =  A(t)e*'^^ 

At  this  point  all  the  disturbance  parameters  have  been  rendered  in  their  time 
dependent  forms,  and  we  reflect  on  the  goal  of  this  analysis.  The  modal  solution  for  the 
perturbation  provides  time  evolution  information  and  leads  toward  stability 
considerations.  To  comment  on  the  stability  of  the  flow  near  the  vortex  sheet  we  need  to 
solve  for  at  least  one  of  the  time  dependent  coefficients  A(r) ,  5j(0 ,  or  82(0 . 
Substituting  the  modal  product  solutions,  (1.5)  and  (1.6),  into  (1.1)  and  rearranging  the 
terms  yields 

(1.7)  = 

and  a  similar  substitution  into  (1.2)  yields 

(1.8)  dfA  =  --^UgikA  +  kB2a!^~^ 

These  relations  give  the  time  rate  of  change  for  the  perturbation  magnitude  and  are 
summed  momentarily.  To  build  a  relationship  between  the  disturbance  potentials,  the 
next  substitution  inserts  the  modal  solutions  into  (1.4)  yielding 

(1 .9)  -  d,B2a'^  -^ikBya'^  -  ^ikB2a'^  =  0 

Pursuing  the  representation  for  A(t) ,  sum  the  time  partials  of  (1.7)  and  (1.8)  giving 
A  =  -^UgikdfA  -  --^UgikdfA  +  ^3,52^*”^ 
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Replacing  the  d,A  terms  in  the  above  sum  with  (1.7)  and  (1.8)  yields 


(1.10) 


+kd,By 


Now  multiply  ka  ’  on  both  sides  of  (1.9)  and  add  the  result  to  (1.10)  removing  all  the 
terms  containing  Bj ’s.  The  resulting  second  order  ODE 


has  solutions  which  are  proportional  to  exponentials 
Ait)oce‘^‘  4-  (T  =  ±^ 

This  result  shows  any  non  constant  perturbation  of  an  infinitesimally  thin  vortex  sheet 
excites  an  instability.  Using  (1.7)  and  (1.8)  we  see  the  stability  of  the  disturbance 
potentials  matches  the  stability  of  the  perturbation.  In  practice,  since  the  vortex  sheet  at 
the  shear  layer  is  not  infinitesimally  thin  and  curvature  cannot  be  ignored  the  instability 
arises  only  for  flows  of  sufficient  Reynolds  number.  Drazin  &  Reid  [12]  address  the 
impacts  of  shear  layer  thickness  and  curvature  on  the  initial  flow  instability.  A  nontrivial 
shear  layer  thickness  directly  affects  the  high  Fourier  modes.  The  stability  of  these 
modes  is  maintained  through  the  viscous  dissipation  within  the  layer.  The  curvature 
impacts  the  low  Fourier  modes.  A  tight  curvature  negates  the  formation  of  low  modes 
since  the  wavelength  cannot  fit  along  the  length  of  the  shear  layer.  A  critical  Reynolds 
number  exists  for  a  given  set  of  flow  parameters  such  that  the  flow  is  stable  when  the 
Reynolds  number  is  below  the  critical  value  and  unstable  otherwise.  After  the  onset  of 
the  instability,  the  flow  equalizes  forming  an  n  vortex  braid  where  n  is  the  Fourier  mode 
with  the  lowest  critical  Reynolds  number.  The  homogeneous  flow  stability  results  in 
section  4  show  the  balance  between  curvature  and  shear  layer  thickness. 

1.2  Problem  Formulation 


The  mathematical  model  for  Newtonian  fluid  flow  has  existed  from  more  than 
150  years.  The  equations  of  motion  describing  physical  fluid  flow  for  this  fluid  are: 
continuity  (mass  conservation),  momentum  conservation,  and  balance  of  energy.  The 
complexity  of  these  equations  necessitates  the  application  of  relevant,  simplifying 
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assumptions  for  the  flow.  The  principal  assumptions  used  for  this  work  are:  2- 
dimentional  flow  in  a  closed  cylindrical  disk,  incompressible/homogeneous  fluid,  low 
velocity  flow,  and  Ekman  forcing  resulting  from  top  and  bottom  friction.  Other 
assumptions  are  discussed  as  they  are  applied  throughout  the  derivations.  Fluid 
continuity  and  momentum  balance  adequately  model  the  flow  field  and  provide  a  set  of 
equations  that  may  be  analyzed  or  solved  numerically.  Since  thermal  excitation  of  the 
flow  is  not  considered  and  buoyancy  effects  are  neglected,  the  energy  balance  does  not 
provide  additional  dynamic  information.  The  time  evolution  of  energy,  enstrophy,  and 
angular  momentum  is  used  to  check  the  numerical  accuracy  of  the  simulation.  Also,  the 
circulation  is  enforced  via  incorporation  into  the  no-slip  boundary  conditions. 

1.2.1  Equations  of  motion 

First  consider  the  continuity  equation.  Continuity  of  a  fluid  is  identical  to 
considering  the  conservation  of  fluid  mass,  and  fluid  continuity  is  represented  as 

(1.11)  $  +  yt>V.«=0 

where  p  is  the  fluid  density  and  u  =  -t-  is  the  fluid  velocity.  The 

continuity  equation  is  transformed  by  the  assumptions  of  fluid  homogeneity,  ^  =  0 ,  and 

incompressibility,  V  •  m  =  0 .  Once  transformed,  (1.11)  simplifies  to  the  incompressibility 
condition 

(1.12)  V-M=0 

i.e.  the  velocity  field  is  solenoidal.  Note:  the  flow  is  assumed  to  take  place  in  only  2- 
dimensions  thus  =  0  leaving  u  =  =l^u^,Ug'). 

The  next  equation  under  consideration  is  the  most  meaningful  for  this  work.  The 
momentum  conservation  equation  (hereafter  referred  to  as  the  Navier-Stokes  equation) 
results  from  Newton’s  law  of  motion  written  for  a  fluid  continuum 

(1.13)  •^  =  — ^Vp  +  V^  +  vV^m-I-jV(V-m) 

where  v  is  the  kinematic  viscosity,  are  conservative  body  forces  per  unit  density, 
and  u  =  (u^,ug)  is  the  velocity  written  in  vector  notation.  The  2D  assumption  neglects 

the  possibility  of  recirculation  in  the  vertical  direction  implying  =  0 .  Considering 
the  terms  on  the  right  hand  side  of  (1.13),  the  first  term  represents  flow  acceleration  due 
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to  a  pressure  gradient  while  the  third  term  shows  the  velocity  dissipation  resulting  from 
viscosity.  The  final  term  vanishes  due  to  the  incompressibility  of  the  fluid. 

Additional  terms  appear  in  the  momentum  balance  when  the  reference  frame  is 
changed  to  a  constant  rotation  in  the  plane  of  fluid  flow.  The  frame  change  removes  a 
rigid  rotation  from  the  flow  and  produces  a  relative  velocity  field 

(1.14)  UQ=u-Qxr 

where  Uq  is  the  velocity  in  the  rotating  frame,  u  is  the  velocity  in  the  fixed  frame, 

Q  =  Qe^  giving  Q  as  the  constant  angular  velocity  of  the  rotating  frame,  and  r  =  re^ 
with  r  being  the  distance  from  the  origin.  The  new  frame  adds  two  accelerations, 
Coriolis  and  centripetal,  into  the  momentum  balance.  The  Coriolis  acceleration  results 
from  the  gradient  in  the  tangential  velocity  and  has  the  form 

(1.15)  2Qxuq  =  +  2Q,^r^e^ 

where  is  the  radial  component  of  the  velocity  in  the  rotating  frame.  The  centripetal 
acceleration  acts  normal  to  the  tangential  velocity  as 

(1.16)  Qx{Qxr)=-Q?re^ 

The  centripetal  acceleration  can  be  written  as  the  gradient  of  a  potential  and  included 
with  the  body  force  term.  Using  vector  identity  A.l  with  a -  b -  Cl  and  c-f 

Clx{pLXf^  =  {pt^^Cl  —  ^f 

the  centripetal  acceleration  equates  to  a  gradient 

(1.17)  Clx{Clxf)  =  -V(l)^ 

where  the  potential  function,  ,  is  given  by 

(1.18)  (Pc=lpxrf 

The  gradient  is  appropriate  in  this  case  since  Q  is  constant  implying  =  d^Q,  =  0 .  As 

a  result  the  gradient  of  the  centripetal  potential  gives 

=-V^jpxr^^=^-Q,^r-Qr^ =-Q^f 

This  leaves  the  Coriolis  acceleration  which  is  discussed  in  more  detail  during  the 
development  of  the  vorticity  form  of  the  Navier-Stokes  equation.  Including  the  above 
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its  magnitude  is  small  when  compared  with  the  primary  2-dimensional  flow.  The  Navier- 
Stokes  equation,  (1.19),  becomes 

(1.21)  d,u  +u  -  Vu  +  2Qxu  =  — kVp  + vV^M  +  Ajd\u 

The  A^d\u  term  results  from  the  vertical  transport  of  the  fluid  due  to  the  Ekman  friction, 

Pedlosky  [24].  Combining  the  2-dimensional  flow  with  the  closed  cylinder  on  top  and 
bottom,  we  assume  the  velocity  profile  is  parabolic  in  the  azimuthal  direction.  Under  this 
assumption  the  vertical  diffusion  is  given  within  the  plane  of  interest  as  a  difference  in 

velocity,  \d\u=^{u*  —U^.  A  complete  discussion  of  this  term  is  provided  in  section 

1.2.2.  Applying  this  representation  for  Ekman  forcing  into  (1.21)  yields  the  final  2D 
approximation  of  the  Navier-Stokes  equation  in  velocity  form 

(1.22)  djM  +  M  •  Vm -(•  2i2  xm  =  Vp -t- -l- —u^ 

where  u  =  (u^,ug'j  is  the  2D  velocity  vector,  p  is  the  fluid  density,  p  is  the  pressure,  v 

is  the  fluid  viscosity,  h  is  the  height  of  the  cylindrical  disk,  and  m*  is  the  forcing  vector. 

With  the  continuity  and  Navier-Stokes  equations  available  in  velocity  form  as 
(1.12)  and  (1.22),  respectively,  we  continue  by  converting  the  Navier-Stokes  equation 
into  the  more  useable  vorticity-stream  function  form.  To  this  end  recall  the  definition  of 
vorticity  as 

(1.23)  cd  =  VxU-a)e^ 

and  the  stream  function  is  defined  by  its  relationship  to  velocity 

(1.24)  u  =  ViffXe^ 

Using  these  definitions  and  the  following  identities,  the  velocity  form  of  the  Navier- 
Stokes  equation  is  transformed  to  the  vorticity-stream  function  form.  Using  vector 
identity  A.2  let  a=b  =  u  yielding 
m-Vm  +3xu 

Using  vector  identity  A.3  let  a -u  yielding 

Vx<j) 

Begin  the  transition  by  substituting  the  above  identities  into  (1.22)  giving 
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(1.25)  d,u  +  (fi)  +  2Q)x  H  =  --I  V/7  - 1  Vm^  -  vV  xft)  +  ^{u  -  u) 

In  order  to  convert  (1.25)  entirely  to  vorticity  and  stream  function,  apply  the  curl  operator 
to  both  sides  of  the  equations.  Applying  the  curl  and  using  the  definitions  (1.23)  and 
(1.24)  gives 

(1.26)  d,a  +  V x((ft)  +  2Q)x (V ^^^x e =  -tA7 x (V x (d)+^{w  - m) 

Notice  the  gradient  fields  vanish  under  the  curl.  One  must  be  careful  since  the  pressure 
has  been  removed  from  the  equations  via  a  mathematical  procedure.  The  physical 
existence  of  the  pressure  cannot  be  ignored.  For  our  problem  with  a  completely  enclosed 
fluid  any  force  due  to  pressure  is  balanced  through  an  equal  force  at  the  boundary. 
However,  analysis  incorporating  a  free  surface  must  adapt  boundary  conditions  to 
resolve  the  pressure  differential  at  the  free  surface.  A  detailed  discussion  of  the  pressure 
calculation  is  beyond  the  scope  of  this  work.  Continuing  the  derivation  using  vector 
identity  A.3  with  5  =  V  x  m  =  d)  yields 

V  X  (V  X  d>)  =  -  V^d) 

while  vector  identity  A.4  with  a  —  a>  and  b  =  Vi//  yields 

V  X  (<yV  X  V 

Applying  these  identities  to  (1.25)  gives 

(1 .27)  +  V  (d>  +  2Q)x  V  +  f{cd*-  (o) 

Now  the  Coriolis  acceleration  is  removed  when  considering  the  vorticity  formulation 
because  Q  is  a  constant  implying  VQ  =  0 .  The  final  form  of  the  Navier-Stokes  equation 
in  vorticity-stream  function  form 

(1.28)  ^,d)  +  7(©,^^r)  =  tA72d)  +  ^(ft)*-ft)) 

where  O)  is  the  magnitude  of  vorticity  vector,  ij/  is  the  stream  function,  and  is 

the  Jacobian,  . 

From  the  standpoint  of  numerical  solvability,  (1.28)  gives  a  Helmholtz  equation, 
and  taking  the  curl  of  (1.24)  produces  a  Poisson  equation.  These  equations 

(1.29)  =  -6))-J{a),\p) 
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(1.30)  V^f  =  -G) 

constitute  the  vorticity-stream  function  formulation  for  this  problem  and  are  the 
mathematical  model  including  all  our  assumptions.  The  method  for  solving  the  system 
(1.29)/(1.30)  is  discussed  in  section  2. 

1.2.2  Ekman  friction 

The  Ekman  friction  on  the  top  and  bottom  of  the  disk  induces  vertical  motion  in 
the  fluid  through  a  phenomenon  known  as  the  Ekman  spiral.  This  motion  occurs  within 
the  Stewartson  layers  and  shown  in  figure  1.3.  The  magnitude  of  the  Ekman  spiral 
velocity  in  the  plane  is  small  when  compared  to  the  primary  rotation.  Thus,  the 

component  of  Ekman  velocity  in  the  plane  is  neglected  while  the  vertical  component  is 
retained.  To  incorporate  the  fluid  pump  created  by  the  Ekman  friction,  the  Ekman 
number  is  defined 


where  v  is  the  viscosity,  Q  is  the  mean  rotation  of  the  system,  and  h  is  the  height  of  the 
fluid  in  the  disk.  The  vertical  velocity  induced  by  the  Ekman  pumping  can  be  written  in 
terms  of  the  planar  vorticity,  and  Friih  [14]  relates  the  vertical  velocity  to  vorticity  using 

where  Ro  =  is  the  Rossby  number,  U  is  the  horizontal  velocity  scale,  and  L  is  the 
Ekman  layer  thickness.  The  Ekman  layer  thickness  is  given  as 
L  =  E^^h 

per  the  Stewartson  layer  thickness  of  vertical  transport.  Similarly,  Pedlosky  [24]  solves 
the  velocity  diffusion  equation  finding 

(1.31)  \dla=^{a*  -a) 

where  \  =  h^E  for  a  geometry  with  a  closed  boundary  on  all  sides.  Both  Friih  [14]  and 

Pedlosky  [24]  performed  asymptotic  expansions  on  the  flow  in  terms  of  the  Rossby 
number.  Combining  these  results  shows  the  vertical  velocity  and  vertical  diffusion  are 
related  through  the  vorticity.  Continuing  with  the  Pedlosky  approach,  assume  the  fluid 
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has  a  parabolic  profile  for  the  velocity  in  the  d-z  plane.  The  horizontal  velocity  scale  is 
given  by 


Adapting  the  horizontal  velocity  scale  into  the  Rossby  number  modifies  the  coefficient  in 
(1.31)  producing 

■Je  _ 

JJ— 

2Q.L  Sa  y  V 

Our  final  assumption  requires  the  disk  height  to  be  close  in  magnitude  to  the  Ekman  layer 
thickness.  This  condition  minimizes  the  influence  of  the  centrifugal  and  Coriolis  forces 
on  the  instabilities.  Thus  entering  h~L  and  the  normalized  distance  to  the  shear  layer, 
a=j  gives  the  form  of  the  Ekman  pumping 

This  term  represents  the  vertical  fluid  motion  as  a  function  of  the  planar  velocity  and 
matches  our  specific  disk  geometry. 

1.2.3  Energy  evolution 

The  equation  for  energy  balance  is  not  used  to  solve  for  vorticity  or  stream 
function,  but  rather  to  provides  a  verification  tool  for  the  numerical  solution.  Thermal 
considerations  are  not  included  with  the  energy  evolution  equation  since  no  external 
thermal  gradients  are  applied  to  the  system  and  the  buoyancy  effects  are  neglected. 
Without  the  thermal  components  the  energy  equation  shows  a  loss  as  viscous  dissipation 
effects  are  converted  to  internal  energy.  The  energy  evolution  equation  is  derived  from 
the  momentum  balance.  Consider  (1.22)  and  dot  the  equation  with  u .  The  result 

•jd, (m •«)+!<•  (m-Vm)  +  m- (2Q X m )  =  -Tn •  Vp  +  VM •  V^M  +^u-(u*  -ii'j 

can  be  integrated  over  the  circular  domain  of  radius  one.  Note:  since  u  lies  in  the  planar 
domain  and  Q  and  co  are  normal  to  the  plane,  the  ensuing  dot  products  are  zero; 

M-i2  =  0,  M-d)  =  0,  m-(<2)xm)  =  0,  and  In  order  to  simplify  the  energy 

balance,  change  to  an  integral  form  by  multiplying  by  a  differential  area  and  integrating 
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(1.32) 


dfjjju  -udS  +  JJm  -(m  ■S/u)dS  = 

-jlju-  WpdS  +  vJJ  M  •  V^udS  +  JJ  u-{u*-  u^dS 

There  is  no  potential  energy  with  this  geometry  and  the  kinetic  energy  is 
T  =  jjju-udS 

Thus  the  time  evolution  of  energy  is  .  To  simplify  (1.32)  consider  each  term 
individually.  The  second  term  on  the  left  hand  side  can  be  rewritten  using  vector  identity 
A.2  letting  a=b=u  and  then  dotting  u  into  both  sides  gives 
M  •  V  (m  •  M  )  =  2m  •  (m  •  Vm  )  +  2m  • 

Using  this  result  the  second  term  is  rewritten  as 

(1.33)  JjM-(M-VM)/5=J'J‘iM-V(M-My5 

Continuing  use  vector  identity  A.5  with  a  =  u-u  and  b=u  yielding 
V  -((m +U-V  {U-u^ 

Applying  this  result  to  (1.33)  gives 

(1.34)  JJm  •  (il •  Vm)/5  =  JJ^V •  ((m  •  u)u)dS 

To  complete  the  final  step  for  this  term,  apply  the  2D  case  of  the  divergence  theorem 
(Gauss’  theorem)  given  here 

JJ  V  •  ddS  =  |J|a  •  fids 

and  (1.34)  becomes 

(1.35)  JJm -(m  •VM)d5  =  |J](m -11)11 -11^5 

where  n  is  the  unit  vector  normal  to  the  boundary  of  the  fluid.  Moving  to  the  right  hand 
side  of  (1.32),  the  first  term  can  be  rewritten  using  vector  identity  A.5.  Let  a  =  p  and 

b  =  u  yielding 

v-(/7m)=  py^ +M -Vp 
and  again  apply  Gauss’  theorem  to  obtain 

(1.36)  JJ  -^M  -  Vpd5  =  JJ  -^V  -  (^pu)dS  =  •  nds 
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Since  our  simulation  is  conducted  in  a  closed  cylinder,  the  unit  normal  vector  is 
everywhere  normal  to  the  surface  of  the  cylinder.  Also,  the  velocity  normal  to  the 
surface  is  necessarily  zero  because  there  is  no  flux  through  the  wall  i.e.  no  fluid  passes 
through  the  cylinder  walls.  Since  velocity  normal  to  the  wall  is  zero,  the  dot  product  with 
the  unit  normal  vector  is  also  zero,  m  •  n  =  0 .  Hence,  both  terms  (1.35)  and  (1.36) 
contribute  nothing  to  the  energy  evolution.  Rewriting  the  energy  equation  shows  that  the 
energy  evolution  results  exclusively  from  the  viscous  effects 

(1.37)  + S’ jj“  ■  (“* 

At  this  point  the  integrands  are  simplified,  and  the  energy  evolution  is  written  for  our 
specific  numerical  work.  Using  vector  identity  A.6,  let  a  =  V  and  b  =  u  then  dot  the 
velocity  into  both  sides  of  the  identity  for 

M  •  -  M  •  =  -M  •[Vx(Vxm)J  =  -m  •  (Vxd>) 

Continuing  use  vector  identity  A.7  with  a  =  m  and  b  =cd  giving 
-M-(Vxd))  =  V-(MXft))-G)(VXd))=  V-(MXd))-fi)-d) 

Combining  the  above  equations  and  employing  Gauss’  divergence  theorem  leads  to 

(1.38)  JjM-V^MrfS=JJV-(MX<2>)/S-W=[[](MXd))-nd5-W 

where  n  is  the  unit  vector  normal  to  the  boundary  and  1^  =  J|  co^dS  is  the  enstrophy. 

Now  the  integrand  of  the  line  integral  is  rewritten 
(m  X  d))  •  n  =  (ugCOe^  -  Uj.coe0  UgCO 

and  is  evaluated  at  the  boundary.  Incorporating  (1.38)  and  expanding  the  second  term  on 
the  right  hand  side  of  (1.37)  produces  an  energy  evolution  with  the  form 

(1.39)  f  =  ‘10-VW  +  ffJ^’u-u'derdr-f-T 

Here  we  note  the  energy  evolution  is  a  function  of  the  enstrophy  and  the  kinetic  energy. 
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1.2.4  Enstrophy  evolution 

While  not  one  of  the  equations  of  motion,  the  enstrophy  evolution  provides  an 
important  diagnostic  test  for  the  numerical  simulation.  The  enstrophy  is  the  mean  square 
measure  of  vorticity,  and  its  evolution  can  by  calculated  analytically.  Beginning  with  the 
Navier-Stokes  equation  in  vorticity  form 

d,a)  =  Yx(ux(&  +  2Qyj+W^(o+^(co* -0)) 

multiply  both  sides  by  the  vorticity  giving 

=di-'^x{ux{B  +  2Q^j^  +  V(b-V^(0-¥^(b-{of  -co) 

As  with  the  energy  evolution,  integrate  this  equation  over  the  unit  circle  giving 

i  JJ  d,a)^dS  =  JJft)  •  [v  X  (m  X  (<j>  +  2i2))]c/5  + 
v|J  Q)  •  V^(bdS  +  •  {(b*  -  m)dS 

Enstrophy  is  defined  by  W  =  JJ  aP'dS  so  the  right  hand  side  of  (1.40)  already  represents 
enstrophy  evolution,  ^ . 

For  the  first  step  to  simplify  the  terms  in  (1.40)  recall  the  2D  assumption  for  the 
flow  in  the  disk.  The  6)  +  2lQ  vector  is  normal  to  the  2D  surface  of  flow, 

6)  +  2Q  -(a)  +  2Q,)e^ .  As  a  result  any  dot  product  between  the  vorticity  vector  and  a 

vector  contained  in  the  2D  flow  equals  zero  e.g.  m  •  (d)  +  2Q)  =  0  and  (d)  +  2Q)  -  V  =  0 . 
Using  the  2D  assumption  and  vector  identity  A.6  with  a  =  u  and  b=d>  +  2Q  yields 
Vx  Xft>  +  2^2^  =  ~  ^d>  +  ~  m  •  V  ^ci)  +  2^2^+  u 

Multiplying  by  d)  and  recalling  2^2  is  a  constant  gives 

d)-[Vx(MXd)  +  2Q)j  =  -d)-(M-Vd))+d)-M- = Vd?^ 

This  expression  matches  the  integrand  of  the  first  term  on  the  right  hand  side  of  (1.40). 
Rewriting  the  first  term  gives 

(1.41)  JJd)  •  [V  X  (m  xd)  +  2Q)]rf5  =  -4  JJ u  ■  V(6^dS 

A 

Next  apply  vector  identity  A.5  letting  a  =  co  and  b=u 
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V  •  (co^u^  =  0)^  +  uVa)^ 

and  Gauss’  theorem  to  (1.41)  giving  the  now  familiar  form 

(1 .42)  JJft)  •  [V  X  (s  X d)  +  JJ V  •  =  -j\l\Ct)^u  ■  fids 

where  n  is  the  unit  vector  normal  to  the  boundary.  Again  the  velocity  normal  to  the 
cylinder  walls  is  zero  negating  any  contribution  from  the  above  term. 

As  with  the  energy  evolution,  only  the  terms  resulting  from  viscous  interactions 
contribute  to  the  enstrophy  evolution.  The  resulting  enstrophy  evolution  equation  is 

(1 .43)  ^  =  2v  JJ  6)  •  +  ^jja)-{co-  Q))dS 

Further  adaptation  of  the  enstrophy  evolution  begins  by  simplifying  the  first  term  on  the 
right  hand  side  of  (1.43).  This  is  accomplished  by  considering  Green’s  theorem 

jjaV^bdS-^\^aWb-n-  ffva  ■VbdS 
while  allowing  a  =  b  =  (0.  The  first  term  is  written  as 
v|J  oN^codS  =  v^oNm  ■  Ms  -  vJJ  (VcofdS 

where  Vtw  =  d^(oe^  +ydgCOe0  in  the  2D  circular  geometry  of  the  disk.  Applying  this 

transformation  and  an  expansion  of  the  last  term  on  the  right  hand  side,  equation  (1.43) 
becomes 

(1 .44)  dR  =  2v^aNa)-  Ms  -  2v JJ {V(of  dS+^\^cb-  mdS  - 

The  evaluation  of  the  integrals  for  our  specific  flow  parameters  begins  with  the 
representation  of  the  integrands 

aNco-n  =  a)(d^(joe^+jdga)e0ye^\^^^  =  -aid^o)\^^^  =-jd^a/ 

_  _  4:  ^  /V  + 

0)  0)  =  coe^  -0)  e^-  0)0) 

The  final  form  for  (1.44)  is 
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1.2.5  Angular  momentum  evolution 

The  angular  momentum  is  another  diagnostic  tool  used  to  check  the  accuracy  of 
the  spectral  simulation.  The  following  discussion  is  adapted  from  Coutsias  &  Lynov  [8]. 
Using  the  given  geometry  we  consider  a  2D  domain  in  polar  coordinates  and  enclosed  on 
all  sides  by  a  rigid  boundary. 

The  scalar  angular  momentum  is  defined  as 

A  =  A-e^  =  JJ(MXr)-ndS 

where  f  is  in  the  plane  and  n  =  e^.  Combining  (1.25)  and  the  above  definition  we  see 
the  time  evolution  of  angular  momentum  is 

f  =  JJ(a,iixr).2.<iS  =  X4 

n=l 

where  the  /„  terms  are 

(1.45)  7i  =JJ[(MXG))xf]-e,d5 

(1.46)  h  =-jj[^{W+Py^}^zdS 

(1.47)  73  =-vJ|[(Vxft))xr]-4fi?5=vJJ(v^MXr)-e^dS 

(1.48)  74 

Each  7„  term  is  evaluated  independently.  For  7j  notice  the  combination  of  cross 
products  gives 

(m  X  d))x  r  =  (ii  •  r  )<5) 

Substituting  the  above  relation  and  the  definitions  for  velocity  and  vorticity  into  (1.45) 
yields 

(1.49)  7i 

The  integrand  is  expanded  to  the  form 

(1.50)  )~de^r )]  =  +  ru^Ug  -  ru^dgU^ 

At  this  point  recall  the  flow  is  divergence  free.  Expanding  the  divergence  of  the  velocity 
while  multiplying  r^Ug  to  both  side  of  the  equation  gives 
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r^Ug  (V  •  m)  =  0  =  rUgU^  +  r^Ugd^U^  +  rUgdgUg 

Adding  the  zero  divergence  into  (1.50)  and  rearranging  the  terms  gives 

[■7(^r  )] “  r^u^d^Ug  +  r^Ugd^u^  +  Iru^Ug  +  rugdgUg  - ru^dgU^ 

=  d,  (r\ug ')  +  dg(rul)-^dg  ) 

Rewriting  the  integral  (1.49)  with  the  above  partial  derivatives  yields 

^  lo  ~  2  Jo 

The  integral  evaluates  to  zero  because  Ug\^^^  =  0  and  =  0 ,  and  the  other  two 
integrals  evaluate  to  zero  because  the  velocity  is  periodic  in  the  6  component.  For  the 
second  term  apply  the  vector  identity  A.4  with  a=^u^  +  p  and  b  =  r  producing 

/2  =  -  JJ  ^  V  (4  +  p  )  X  r  ]  • 

=  V  X  (4m^  +  p)r  •  e^dS  +  JJ  +  p)(^^  ■  e^dS 
Continue  by  applying  Stokes’  theorem 

JJ(Vx  5)- ndS  =  ijd  • 

and  the  second  term  equates  to  zero 

4  =  -JJ^Vx  (4^2  +  p)r]-  =  |](4m^  +  p)r  •  ^  =  0 

since  the  radial  vector  is  normal  to  the  boundary,  7  Lds .  Moving  to  the  third  term  begin 
with  vector  identity  A.3.  Allowing  d-u  while  including  the  definition  for  vorticity 
gives 

V^M  =  -  Vx  (V  xii)  = -V  xd) 

Incorporating  the  curl  of  the  vorticity  and  the  above  identity  with  (1.47)  yields 

73  =  V JJ  X  r  ^  •  e^dS  =  “V' JJ  ~  J^  Jq  f^^r^drdO 

Evaluating  the  integral  gives  the  following  formula 

Ij  =  -IttvcOq  (1) + 4;rvMo  (l) 

where  cOq  and  Uq  are  the  zero  Fourier  mode  terms  for  vorticity  and  velocity,  respectively. 
The  non-zero  Fourier  modes  integrate  to  zero  over  the  range  0 : 0  — >  2;r .  The  final  term 
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is  expanded  for  the  integral.  We  can  reduce  the  integrand  of  (1.48)  since  both  r  and  u 
reside  in  the  plane  giving 

Recall  the  forcing,  ,  is  independent  of  6  and  velocity  is  written  through  the  stream 
function 

Ug  =  d,l/f 

Again  due  to  periodicity  in  the  azimuthal  direction  all  the  non-zero  Fourier  modes  vanish. 
Evaluating  the  final  integral  gives 

Combining  all  of  the  terms  gives  the  angular  momentum  evolution  as 
f  -  =  -IttvcOq  (1) + 4;rVHo  (0 "  ^  [^e  (^)  "  ^rVo  ir)]dr 

Due  to  the  viscosity  the  angular  momentum  is  not  a  conserved  quantity. 

1.2.6  Circulation  evolution 

The  circulation  is  the  measure  of  the  vorticity  flux  through  a  surface  formed  by  a 
closed  curve.  Using  Green’s  theorem  the  vorticity  flux  equates  to  the  component  of 
velocity  tangent  to  the  closed  curve  along  the  closed  curve 

C  =  IJ  oydS  =  ^U‘dl 

For  the  no-slip  boundary  conditions  the  velocity  at  the  wall  is  zero  in  all  directions 
indicating  the  circulation  is  also  zero.  To  fully  investigate  the  circulation  and  its  time 
evolution  begin  with  the  velocity-vorticity  equation 

d,M  =  Vx(uxco)+vV^co+^(a>* 

Integrating  the  vorticity  equation  over  the  unit  circle  immediately  yields  the  circulation 
evolution 

(1.51)  ^  =  jjvx{uxw)-e^dS+  v JJ V^codS  + JJ {co*  - cd)dS 

Employ  Stokes’  theorem 

J|(V  X  a)  •  fids  =  ||]a  •  ds 

with  a=uXQ}  to  rewrite  the  first  term  in  ( 1 .5 1)  as 
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J|Vx(MX«;)-e^c/5  =  ^(uxa>)-dl 

This  term  vanishes  for  either  the  no-slip  or  free-slip  boundary  conditions  since  the  fluid  is 
contained  in  a  rigid  disk.  With  the  no-slip  boundary  conditions,  the  velocity  is  identically 
zero  implying  the  line  integral  is  also  zero.  Since  there  is  no  flux  through  the  boundary, 
the  free-slip  conditions  give  a  boundary  velocity  that  is  tangent  to  the  wall,  . 

Applying  the  curl  between  this  velocity  and  the  vorticity,  u^eg  x  a)e^ ,  results  in  a  vector 
directed  along  the  r  axis.  Since  dl  =  dlcg ,  the  dot  product  with  the  curl  of  the  integrand 
and  the  differential  length  equals  zero,  UgO)e^  ■  dleg  =  0 .  Thus  the  term  equates  to  zero 

for  both  sets  of  boundary  conditions.  Now  apply  the  divergence  theorem  to  the  second 
term  on  the  right  hand  side  of  (1.51)  giving 

vJJ  V^a)dS  =  vJJ  V  •  V  cods  =v^S/co-ndl 

Combining  this  new  form  with  a  rewrite  of  the  last  term  on  the  right  hand  side  of  the 
circulation  evolution  results  in  the  form 

(1.52)  f  =  •  ndl  +^J|  a>*dS  ~fC 

Similar  to  the  enstrophy  evolution,  the  time  rate  of  change  for  circulation  depends 
explicitly  on  the  current  value  of  the  circulation.  This  effect  results  directly  from  the 
inclusion  of  Ekman  friction.  Before  solving  this  ODE  we  consider  a  final  transformation 
for  the  specific  values  of  the  integrals  on  the  unit  disk.  For  the  first  term  on  the  right 
hand  side  of  the  circulation  evolution 

Vft>  • «  =  +  ^dffCOeg  )  •  {-e^ )  =  -d,  co\^^^ 

For  the  second,  the  forcing,  d)* ,  is  a  function  of  r  only  and  is  independent  of  6 . 
Combining  these  ideas  into  (1.52)  gives 

(1.53)  f  = 

where  cOq  is  the  constant  term  from  the  Fourier  series  expansion  of  the  vorticity.  All  the 

higher  Fourier  modes  vanish  as  the  integral  for  is  evaluated.  When  Ekman 

friction  is  not  considered  with  the  Navier-Stokes  equation,  the  circulation  evolution  is 
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zero.  To  maintain  consistency  the  value  of  is  set  equal  to  zero  in  the  numerical 

simulation.  The  ODE  provided  by  (1.53)  gives 


showing  the  circulation  decays  exponentially  to  a  steady  value  consistent  with  the 
forcing. 

1.2.7  Dimensional  analysis  and  scaling 

For  ease  of  expanding  the  functions  in  terms  of  Chebyshev  polynomials,  the 
dimensionless  radial  variable  is  scaled  to  [-1,1].  The  actual  radial  distance  exists  on  the 
interval  0<r<R  so  the  dimensionless  distance  on  -1  <  r  <  1  results  from  the  scaling 
r  =  -j.  The  angular  variation  is  dimensionless  and  occurs  on  [0,2;r] .  The  full  list  of 
dimensionless  parameters,  shown  with  a  ~,  is  given  here 


Su  =  aSQ, 


where  Su  is  the  differential  tangential  velocity,  a  is  the  radial  distance  to  the  shear  layer, 
^  is  the  differential  angular  velocity,  and  P  is  the  appropriate  length  scale  for  the 
Reynolds  number.  For  consistency  the  Ekman  forcing  vorticity  is  made  dimensionless 
via 

(1-54) 

to  preserve  the  equality  -0)^  =  ~ •  ^f^er  applying  the  dimensionless 

parameters  to  (1.28),  the  result 

(1.55)  djQ)  + =  +  -<y)j 

can  be  rewritten  using  the  Reynolds  number 
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(1.56)  d-^(b  +  J{cb,xj/)  =  ^^^cb+Y^^*~^) 

The  Reynolds  number  takes  on  the  standard  form  for  fluid  flow  in  a  circular  geometry 
where  the  Kelvin-Helmholtz  instability  is  studied 

(1.57)  Re  =  -^ 

Equation  (1.30)  maintains  identical  form  after  applying  the  dimensionless  parameters 

(1.58)  vV  =  -fi> 

From  this  point  on  the  ~  is  dropped  and  all  variables  are  understood  to  be  dimensionless. 
The  most  important  adjustment  for  future  calculations  occurs  with  the  forcing,  (1.54). 
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2  Numerics 

2.1  Introduction 

Spectral  methods  are  employed  for  numerical  simulation  by  expanding  the 
variables  of  interest  using  complete  sets  of  global,  orthogonal  functions.  The  global 
nature  of  spectral  methods  make  them  ideal  for  solving  partial  differential  equations 
(PDEs)  in  simple  geometries.  Fomberg  [13]  discusses  the  advantages  of  expanding 
variables  using  global  anal)'tic  functions:  1)  errors  typically  decay  at  exponential  rather 
than  pol)momial  rates,  2)  the  dissipative  and  dispersive  errors  are  virtually  nonexistent, 
and  3)  computational  grids  are  typically  courser  allowing  for  faster  computation  and 
lower  memory  requirements.  Some  of  these  advantages  are  lost,  however,  when  spectral 
methods  are  extended  to  complex  domains,  Gottlieb  &  Orzag  [16]. 

Once  an  appropriate  set  of  basis  functions  are  selected,  there  are  three  main 
techniques  used  to  determine  the  expansion  coefficients:  Galerkin  methods,  tau  methods, 
and  pseudospectral  (collocation)  methods.  The  Galerkin  technique  recombines  the  basis 
functions  to  form  a  new  basis  which  automatically  satisfies  the  boundary  conditions.  The 
tau  technique  selects  the  expansion  coefficients  that  satisfy  the  boundary  conditions  while 
minimizing  the  residual.  Our  simulation  uses  the  pseudospectral  technique  where  the 
nonlinear  terms  are  calculated  at  the  grid  points  while  spatial  derivatives  are  computed  in 
the  spectral  domain.  Differential  operators  are  transformed  into  well-conditioned  integral 
operators  with  boundary  conditions  enforced  as  tau  constraints  on  the  expansion 
coefficients. 

Here  we  solve  the  2D  Navier-Stokes  equations  in  a  circular  geometry  referred  to 
as  the  “disk”.  The  disk  is  modeled  with  a  small  height  compared  to  its  radius  resulting  in 
quasi  2-dimensional  (Q2D)  fluid  motion.  Equations  describing  the  flow  are  given  in 
section  1.2.1.  The  basis  functions  used  for  the  disk  are  Fourier  in  the  periodic,  angular 
direction  and  Chebyshev  in  the  bounded,  radial  direction.  The  Fourier  and  Chebyshev 
bases  are  selected  due  to  the  existence  of  fast  transforms  which  are  required  to  evaluate 
nonlinear  terms  in  a  collocation  method.  The  computational  grid  for  the  domain  is 
equally  spaced  in  the  azimuthal  direction  and  uses  Guass-Lobatto  points  in  the  radial 
direction.  The  spatial  variables,  vorticity  and  stream  function,  are  represented  by 


27 


Fourier-Chebyshev  expansions,  and  the  general  expansion  of  these  flow  parameters  has 
the  form 

M 

(2.1)  i(r,0)=  £  -R<r<Ji  0<9<27t 

where  ^  is  either  vorticity, O) ,  or  stream  function,^ ;  {r,9')  is  the  radial  and  angular 
location  in  the  disk,  respectively;  T,„{r)  is  the  m‘^  Chebyshev  polynomial;  is  the 
Fourier-Chebyshev  expansion  coefficient;  y  -I- 1  is  the  number  of  Fourier  modes;  and 

M  + 1  is  the  number  of  Chebyshev  modes. 

The  Fourier  modes  of  the  flow  parameters  uncouple  and  the  resulting  system  of 
uncoupled,  time  dependent  ODEs  is  solved  using  a  third  order,  stiffly-stable  semi- 
implicit  Backward  Difference  Formula  (SBDF).  This  scheme  uses  an  implicit 
formulation  of  the  vorticity  to  solve  the  linear  terms  and  an  explicit  formulation  for  the 
nonlinear  terms 

La)  =  j((o,y/)+  F 

The  SBDF  scheme  is  implemented  with  the  Chebyshev  tau  method  producing  a 
Helmholtz  operator.  Postconditioning  of  the  solution  provides  a  banded  system,  and 
when  combined  with  the  fast  transforms,  yields  a  numerical  solution  in  O  (ATM  log  NM  ) 
operations,  Torres  &  Coutsias  [32]. 

Using  polar  coordinates  on  the  disk  produces  a  point  singularity  at  the  origin,  and 
methods  for  treating  the  singularity  have  been  investigated  by  several  authors.  For  our 
implementation  an  multiplier  is  included  with  the  vorticity  building  partial  regularity 
into  the  flow  field.  A  more  elegant  solution  to  the  point  singularity  involves  an 

orthogonal  basis  behaving  as  O  {r'"  ^  as  r  — >  0 .  This  method  is  treated  by  Matsushima  & 

Marcus  [22]  but  lacks  the  fast  transform  required  for  numerical  simulation.  Several 
boundary  condition  solutions  are  offered  to  treat  the  singularity,  and  these  can  be  found 
in  Torres  &  Coutsias  [32]  and  the  included  references. 

The  code  written  in  Fortran  90  is  paralleled  using  the  Message  Passage  Interface 
(MPI)  parallel  language.  The  computational  limitations  arise  due  to  the  Fast  Fourier 
Transform  (FFT).  During  the  transposition  of  the  matrices  used  in  the  calculations, 
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information  is  passed  between  paralleled  machines  limiting  computational  speed  to  the 
computer  interface  hardware.  A  table  showing  the  computational  efficiency  of  the  code 
can  be  found  in  Torres  &  Coutsias  [32]. 

2.2  Fourier-Chebyshev  expansion 

Using  a  spectral  scheme  within  the  disk  geometry  lends  itself  to  the  Fourier- 
Chebyshev  expansion  of  the  stream  function  and  vorticity.  The  Fourier  expansion  in  the 
azimuthal  direction,  6 ,  matches  the  natural  periodicity  resulting  from  polar  coordinates. 
The  Chebyshev  expansion  in  the  radial  direction,  r ,  allows  for  efficient  postconditioning 
of  the  problem  and  fast  transformations  between  point  space  and  mode  space.  The 
resulting  expansions  for  stream  function  and  vorticity  follow  from  (2.1)  and  are  given 
here  considering  a  trigonometric  Fourier  basis 

■f"  M 

(2.2)  =  + 

n=0  m=0 

^  HA 

2  M 

(2.3)  (O{r,0)  =  KnT,„ir)cos  nd  +  0)^,^  (r)  sin  nO 

/i=0m=0 

where  ,  and  are  the  cosine  and  sine  coefficients  for  the  stream 

function  and  vorticity,  respectively,  -R  <r<R,  and  0<6<27r .  Both  Fourier  and 
Chebyshev  bases  are  orthogonal  and  possess  the  Fast  Fourier  Transform  and  Fast  Cosine 
Transforms  (FFT/FCT),  respectively.  The  fast  transforms  offer  O  (AM  log  NM  ) 

operations  to  move  between  point  space  and  mode  space  making  them  ideal  for 
numerical  simulation.  Once  the  problem  is  transformed  to  mode  space  the  nonlinear 
terms  are  represented  by  convolutions.  The  number  of  computations  required  to  resolve 

the  convolutions  increases  quadratically  with  the  number  of  modes,  (AM  )  ,  destroying 
the  AM  log  AM  pace  of  the  computations.  Therefore,  the  FFT/FCT  are  required  for 
calculating  the  nonlinear  terms  by  transforming  them  into  point  space  where  they  are 
multiplied  directly  in  AM  computations. 

The  grid  applied  to  the  computational  domain  expands  the  Itt  radians  equally 
and  the  individual  grid  points  are 
0^  for  0<A:<A 
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The  radial  Chebyshev  polynomials  are  expanded  on  -R<r<R  instead  of  0  <  i?  using 
the  Gauss-Lobatto  quadrature  points 
^  IcTT 

4=/?cos  —  forO<^<M 

This  choice  of  grid  improves  the  condition  numbers  for  the  spectral  matrices,  Torres  & 
Coutsias  [32].  It  must  be  noted  that  an  expansion  of  0<6<2n  and  -R<r<R 
produces  a  grid  that  covers  the  disk  twice.  The  solution  on  0<6<7r,  r:R-^-R 
matches  exactly  with  the  solution  on  ;z  <  6  <  2;z ,  r:-R-^R.  Since  the  Fourier  modes 
are  uncoupled,  an  important  result  of  this  double  coverage  appears  in  the  spectral 
matrices.  Consider  a  function  on  the  disk 
A(r,6)  =  f(r)cosn0 

where  /(r)  is  expanded  in  Chebyshev  polynomials.  Since  the  gird  on  the  disk  overlaps 
itself  exactly,  the  correlation  of  grid  points  gives  the  following  equivalence 
A{r,6)  =  f{r)cosn0  =  f(—r)cos(0  +  7r)n  =  A(—r,0  +  tt) 

Continuing  f{-r)cos{0+7i:)n  =  f(-r)cosn0(-l)"  implies 

and  the  terms  in  the  Chebyshev  expansions  have  the  same  even/odd  parity  as  the  Fourier 
mode  number,  n.  Furthermore,  f(r)smn0  yields  an  identical  parity.  The  vorticity  and 
stream  function  are  represented  by  a  matrix  containing  the  Fourier-Chebyshev 
coefficients.  The  Chebyshev  polynomials  are  identified  by  the  rows  in  order  0  — >  M  , 
and  the  Fourier  trigonometric  functions  appear  as  columns  paired  cosine/sine  for  modes 
0  — >  Y .  The  parity  existing  between  the  Chebyshev  and  Fourier  mode  numbers  implies 
half  of  the  coefficient  matrix  is  always  zero,  and  the  form  of  the  matrix  is  given  here 


He 

0 

0 

0 

* 

•  •  •  •  •  *  ^ 

* 

0 

0 

* 

* 

0 

0 

.  0 

0 

* 

0 

0 

0 

* 

•  •  •  •  •  •  ^ 

0 

0 

* 

* 

0 

0 

.  0 

0 

0 

0 

0 

* 

.  ^ 

* 

0 

0 

* 

* 

0 

0 

.  0 

0 

* 

0 

0 

0 

* 

.  * 

30 


where  *  represents  a  nonzero  coefficient.  This  guaranteed  structure  for  the  vorticity  and 
stream  function  matrices  allows  a  significant  reduction  in  the  number  of  computations 
since  the  nontrivial  information  can  be  packed  into  a  matrix  of  only  half  size. 

2.3  Spatial  differentiation  in  the  spectral  scheme 

Herein  lies  the  power  behind  spectral  methods.  The  system  of  equations  to  be 
solved  is  dominated  by  radial  derivatives.  The  Chebyshev  polynomials  are  differentiated 
exactly  so  that  the  differentiation  error  derives  from  truncation  of  the  expansion.  Due  to 
the  convergence  of  the  Chebyshev  expansion,  the  error  inherent  in  the  Chebyshev 
expansion  decreases  exponentially  as  more  polynomials  are  retained.  Since  the  same  is 
true  of  the  Fourier  expansion,  computation  of  spatial  derivatives  is  exponentially  accurate 
when  integration  conditioning  is  used.  Two  spectral  matrices  are  constructed 
representing  linear  multiplication  and  differentiation  operators  within  the  Chebyshev 
basis,  see  appendix  D.  The  linear  multiplication  matrix,  R ,  and  differentiation  matrix, 

D ,  apply  to  a  column  vector  consisting  of  Chebyshev  coefficients. 

To  address  the  full  problem,  the  coefficient  matrix  for  vorticity  is  rewritten  as  a 
single  column  vector.  Each  column  in  the  vorticity  matrix  contains  the  Chebyshev 
expansion  coefficients  for  a  specific  Fourier  mode.  These  columns  are  stacked  on  top  of 
one  another  to  form  the  following  vectors  for  vorticity 


cooir) 

d)[(r) 

G)i(r) 

&(r,e)= 

ft)f(r) 

d>i(r) 

2 


nir) 

fW 

mr) 

and  for  stream  function;  \p^{r,9)=  ^ 

0 

Wlir) 

flir) 

2 


where  G)J(r)  and  d)J(r)  are  the  vorticity  Chebyshev  coefficient  vectors  for  the 


Fourier  mode  cosine  and  sine,  respectively,  and  similarly  and  ^j(r)  are  the 

stream  function  Chebyshev  coefficient  vectors.  This  formulation  of  the  problem  is 
hereafter  referred  to  as  the  full  vector  form.  Since  multiplication  and  differentiation  by  r 
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impact  only  the  Chebyshev  coefficients,  the  matrices  providing  multiplication  and 
differentiation  by  r  for  the  full  vector  problem  are 

R  =  0R 

2 

D  =  I;v^,  ®D 
2 


where  is  the  +  +  identity  matrix.  In  addition  to  these  basic  operations 


within  the  Fourier-Chebyshev  expansion,  R  =  r  and  D  =  0,. ,  Fourier  modes  require 
differentiation  hy  d .  In  the  full  vector  environment  a  new  matrix  is  defined  for  this 
purpose,  Q  =  .  The  Q  matrix  is  formed  by  considering  the  differentiation  of  sine  and 

cosine  functions.  Here  the  impact  is  strictly  upon  the  Fourier  modes  leaving  the 
Chebyshev  expansions  unchanged.  For  differentiation  by  d 

0 


Q  = 


0 


where  n  - 1  changes  as  the  row  of  1^,^, .  The  complete  derivations  of  R ,  D ,  and  Q  are 

presented  in  appendix  D.  For  the  numerical  simulation,  differentiation  by  6  only 
appears  as  a  second  partial,  8^ .  The  matrix  form  of  this  second  partial  derivative,  .  is 
given  by 


Q  =1. 


-nh.. 


The  block  diagonal  form  of  R ,  D ,  and  Q  allows  uncoupling  of  the  Poisson  and 
Helmholtz  equations  into  the  Chebyshev  polynomials.  Consider  the  Poisson  equation 

and  multiply  by  to  remove  the  polar  singularity.  Rewriting  in  matrix  representation 
gives 

R^D^+RD  +  Q^)vr  =  -R26) 

Similarly  the  Helmholtz  equation 

d,o) - V )co  =  -d^\i/dgCo)+^{(o*  -(o) 

multiplied  by  and  writing  in  matrix  form  gives 
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(2.5) 


R^a,6)  -  V  (r^D^  +1tD  +  Q^—^R^)ci>  = 
R(Ry/*Qm-RQ)*Qi^^-^R^&* 

The  nonlinear  terms,  R  *  Qa)  and  R<j)  *  Q  ,  are  determined  in  point  space  with  the 
results  transformed  into  a  modal  matrix  via  the  FFT/FCT.  Both  (2.4)  and  (2.5)  are 
uncoupled  and  for  each  Fourier  mode  n  the  system  of  equations  has  the  form 

(2.6)  (R^D^  +  RD  -  =  -R^^„ 

(2.7)  R^d,Q)„  -v(r"d"  +RD-n"l-^R")d)„  =RV„(d),(^)-^R2d>,: 
where  and  co„  are  vectors  containing  Chebyshev  coefficients  for  Fourier  mode  n , 
J„(Q),ip)  is  the  Jacobian  previously  computed  in  point  space,  and  the  forcing  ft>*  =  0  for 

n  >  0.  Fast  computations  are  obtained  for  the  Poisson  equation,  (2.6),  since  R  is 
tridiagonal  and  D  is  postconditioned  by  P  resulting  in  a  banded  system.  The  Helmholtz 
equation  (2.7)  is  formed  using  the  SBDF  method  described  in  section  2.5.5. 

2.4  Postconditioning 

The  spectral  matrices  resulting  from  the  Chebyshev  basis  provide  a  problem 
poorly  conditioned  for  numerical  computation.  The  entries  in  the  matrix  grow  with 

the  Chebyshev  mode  number  at  the  rate  .  A  postconditioner  (or  alternatively  a 
preconditioner)  not  only  improves  the  conditioning  of  the  problem  but  produces  a  banded 
system  of  equations  by  transferring  the  differentiation  operator  into  an  equivalent  integral 
operator,  Torres  &  Coutsias  [32].  The  postconditioner,  P ,  is  based  upon  the  pseudo 
inverse  of  D  such  that  the  matrix  product  DP  yields  the  identity  matrix  whose  rows 
have  been  shifted  up  by  one 

'0  1  0  •••  O' 

0  0  1  ■•.  : 

DP=  :  :  ■•.  ■•.  0 

0  0  •••  0  1 

0  0  •••  0  0 

Similarly  D^P^  gives  the  identity  with  the  rows  shifted  up  by  two.  The  stream  function 
is  modified  with  the  postconditioner 
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which  is  efficiently  implemented  in  practice  since  is  pentadiagonal.  Consider  (2.6) 
with  the  postconditioner  applied  to  the  stream  function 

(2.8)  (R^Ip,  +  RI[,|P  -  )#;  =  -R^S„ 

where  is  the  identity  matrix  shift  up  j  rows.  Boundary  conditions  are  the  final 

consideration  required  to  solve  (2.8).  The  r  conditions  are  implemented  by  shifting  all 
of  the  matrix  entries  down  two  rows  and  inserting  the  boundary  conditions  as  the  first 

two  rows.  The  equation  augmented  with  boundary  conditions  produces  which  is 

easily  converted  back  to  y/^  through  multiplication  by  the  banded  matrix,  . 

For  the  Helmholtz  equation  the  spectral  accuracy  of  the  postconditioned  method 
solution  degrades  as  more  Fourier  modes  are  included.  Alleviating  the  degradation  is 
accomplished  by  smoothing  the  vorticity  during  postconditioning.  Equation  (2.7)  is 
treated  separately  for  n  =  0 ,  n  =  1 ,  and  n>2.  Each  case  is  treated  with  a  different 
postconditioning  of  the  vorticity  aimed  at  smoothing  the  functions  near  the  singularity. 
For  the  case  n  =  0,  transformation  of  the  vorticity  by  a  radial  multiplier  and 

postconditioner  gives  3,.d)o  =  RP^ft)o .  Similarly  for  the  cases  where  n  =  l  and  n>2  the 
postconditioning  comes  in  the  forms  d),  =  RP^d)i  and  d)„  =  R^P^d)° ,  respectively. 
Applying  the  radial  multiplier,  R  or  R^ ,  builds  partial  regularity  into  the  solution,  but 

the  multiplication  is  valid  since  each  Fourier  mode,  n ,  decays  like  r"  as  r  — >  0 . 
Smoothing  the  function  successfully  provides  degradation  abatement  and  is  fully 
analyzed  in  Torres  &  Coutsias  [32].  Applying  the  most  frequent  case,  n  >  2 ,  of  the 
postconditioning  to  the  Helmholtz  equation  giving 

R^P^3,ft);  -v(r^D^RV  +RDRV  -/2^r¥^  +^R^P^)d),;  = 

R^J„icd,w)-fR^cdl 

This  equation  simplifies  through  the  relationship  DR  =  RD  + 1 .  After  simplifying  and 
removing  the  common  terms  we  have 

(2.9)  R^p^a,s;--'[R%2]  +  5Ri[,]P+(4-n")p"+|-i]s: 
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The  time  differentiation  is  expanded  via  the  SBDF  giving  an  equation  for  the 
postconditioned  vorticity  at  the  next  time  level,  .  As  with  the  Poisson  equation,  the 
entries  in  the  matrices  are  shifted  down  two  rows  and  the  boundary  conditions  are 
inserted  at  the  top.  After  calculating  the  vorticity  is  quickly  reclaimed  through  the 

/  _L  t  O  O  /  _1_  1 

multiplication  by  banded  matrices,  =  R  P  d)°’  .  Fourier  modes  0  and  1  follow  a 

similar  development  but  require  special  treatment.  Fourier  mode  1  takes  the 
transformation  =  RP^d)j while  zero  Fourier  mode  is  =  RP^d)o  . 

2.5  Computational  algorithm 

This  section  outlines  the  primary  steps  taken  in  the  code  and  the  order  in  which 
they  occur.  As  a  summary  tool,  each  step  identifies  subsections  relating  to  that  step. 

1)  The  random  noise  is  generated  at  the  grid  points  see  2.5.1 

2)  The  forcing  of  the  shear  layer  combines  with  the  random  noise  to  make  the 
initial  field  see  2.5.2 

3)  The  initial  field  is  projected  onto  the  vorticity  solvability  constraints  with  a 
high  order  adjustment  see  2.5.3 

4)  The  vorticity  is  transferred  to  the  start-up  subroutine 

5)  The  equations  are  solved  in  the  start-up  routine  as  explained  in  steps  6)  - 
10).  Vorticity  and  the  Jacobian  are  transferred  back  to  the  main  code 
see  2.5.4 

6)  The  current  vorticity  drives  the  Poisson  equation  which  is  solved  for  the 
stream  function  see  2.5.3  and  2.5.5 

7)  The  nonlinear  terms  appearing  in  the  Jacobian  are  determined  using  the 
current  stream  function  and  vorticity 

8)  The  forcing  and  Jacobian  are  combined  to  form  the  explicit  portion  of  the 
SBDF  see  2.5.2 

9)  The  Helmholtz  equation  is  formed  by  applying  the  SBDF  scheme  to  the 
vorticity  and  the  explicit  terms  found  in  step  8)  see  2.5.3  and  2.5.5 

10)  The  Helmholtz  equation  is  solved  for  the  new  vorticity 

1 1)  The  diagnostics  are  occasionally  invoked  and  the  relative  errors  of  energy, 
enstrophy,  and  angular  momentum  evolution  are  calculated  see  2.5.6 
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12)  steps  6)  -  1 1)  are  repeated  for  a  specified  number  of  iterations 

2.5.1  Initial  conditions 

The  initial  conditions  used  for  these  computations  are  double  periodic  without 
circulation.  Due  to  the  difficulty  of  producing  a  random  field  on  a  non  uniform  grid,  the 
random  initial  conditions  are  established  on  a  square  where  the  grid  points  forming  the 
square  are  equally  spaced  in  both  directions.  The  periodicity  of  the  initial  data  is 
enforced  in  each  direction  on  the  grid  by  considering  a  random  walk  on  the  unit  circle  as 

where  ^  is  a  uniformly  random  number  between  0  and  1  selected  separately  for  each 
Xj  „  and  y j  „ ,  respectively.  The  noise  enters  as  perturbations  to  the  Fourier  modes.  The 

square  grid  consists  of  four  point  clusters  that  ensure  dual  periodicity 

cos  cos  yj  „  sin  cos  y ^ 
cos  sin  y^„  sin  sin  y 

where  the  value  at  every  point  is  weighted  with  a  decaying  Gaussian.  The  final  form  for 
the  initial  noise  on  the  square  grid  is 
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where  Aj  =  cosxj^ ,  Aj  =  y},!c 

.  _  -2  •  __  .2 

=  e  and  w’y=&  ^  being  the  appropriate  Gaussian  weighting  at  each  location. 

After  the  random  data  are  formed  on  the  square,  the  data  are  interpolated  onto  the 
computational  grid.  The  resulting  interpolation  provides  an  initial  disturbance  without 
bias  from  the  non  uniform  grid  spacing. 
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2.5.2  Forcing  of  the  flow 

The  force  transmitted  to  the  fluid  through  the  Ekman  pumping  acts  as  rigid  body 
rotation  near  the  center  of  the  disk  while  vanishing  at  the  outer  wall.  A  tangent 
hyperbolic  function  models  the  radial  velocity  profile  giving  an  infinitely  differentiable 
flow  near  the  shear  layer.  The  forcing  is  uniform  in  the  azimuthal  direction.  This  forcing 
imparts  a  velocity  only  in  the  angular  direction 

A  plot  of  the  velocity  versus  radius  for  any  angle  shows  the  transition  of  the  velocity  at 
the  shear  layer,  see  figure  2.1.  Since  the  code  solves  the  Navier-Stokes  equations  in 
vorticity-stream  function  form,  the  forcing  is  entered  as  a  vorticity 

a*  =^ri  +  tanh^-Jisech2^1 

A  plot  of  the  vorticity  at  the  shear  layer  is  shown  in  figure  2.2.  With  the  no-slip  boundary 
conditions,  the  initial  vorticity  is  projected  into  the  solution  space.  A  random  vorticity 
field  does  not  necessarily  conform  to  the  solvability  constraints  imposed  on  the  vorticity 
by  the  no-slip  conditions.  The  projection  is  accomplished  by  adding  a  radial  power 
function  to  each  Fourier  mode  in  the  vorticity  expansion 

K{r)  =  0)l{r)+'nnr^ 

where  p  matches  the  parity  of  the  Fourier  mode,  n ,  and  is  selected  as  a  large  number. 
The  coefficient  is  defined  by 

_ 

6„RV'’ 

where  is  the  Chebyshev  expansion  vector  for  .  The  typical  order  of  magnitude 

adjustment  for  projecting  the  vorticity  onto  the  solvability  constraints  is  which 
equals  the  noise  resulting  from  the  FFT/FCT. 
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Velocity  at  the  Shear  Layer 

n  or 

mm 

.6 

\J 

P 

“U.Zj 

/ 

-U.J 

A  nc  _ 

/ 

-U./J 

1  - 

f 

-1 

1  - 

7 

-l.Zj 

1  <  - 

/ 

-l.D 

/ 

-1.  /j 

-2  - 

r 

-Z.Zj  " 

0 

.4  0 

.4  0 

.4  0 

1  ( 

.4  0 

3  i 

I 

.4  0 

ladii 

.5  0 

as  (r 

.5  0 

) 

.5  0 

i  ( 

.5  0 

3  J 

.5  0 

1 

fig  2, 1  Velocity  profile  at  the  shear  layer,  AQ  =  1 ,  a  =  0.5 ,  and  Ar  =  0.05 . 


radius  shown  at  the  shear  layer,  AQ  =  1 ,  a  =  0.5 ,  and  Ar  =  0.05 . 
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2.5.3  Boundary  Conditions 

While  the  code  can  implement  free-slip  and  no-slip  boundary  conditions  (BCs), 
only  the  no-slip  conditions  have  been  used  here.  The  disk  physically  terminates  in  all 
directions  creating  a  solid  boundary  enclosing  the  fluid.  Recall  the  definition  of  velocity 
related  to  the  stream  function 


u  =  -I-  UgCg  =  \de\l/e,  -  d^ifre^ 

Since  the  walls  of  the  enclosure  are  solid,  no  fluid  flux  exists  at  the  boundary  implying 
the  velocity  normal  to  the  wall  is  zero,  =  0 .  Note:  during  implementation  the  outer 
radius  of  the  disk  boundary,  R ,  is  normalized  to  one,  /?  =  1 .  The  no-slip  condition 
implies  a  fixed  velocity  tangent  to  the  wall,  =  U  •  These  condition  produce 

Dirichlet  and  Neumann  constraints  on  the  stream  function.  The  zero  flux  requirement 
gives 


while  the  Neumann  conditions  result  from  the  no-slip  condition 
^r¥n\r=l  =0  for  «  >  0 

=U 


(2.10) 


Combining  the  zero  partials  and  the  potential  nature  of  the  stream  function  gives  the 
Dirichlet  conditions 


(2.11)  ¥„U=0  forn>0 

The  combination  of  Dirichlet  and  Neumann  conditions  on  the  stream  function  means  the 
Poisson  equation  is  overdetermined  and  the  Helmholtz  equation  underdetermined. 

Solving  the  entire  system  requires  transferring  one  set  of  BCs  onto  the  vorticity.  For  our 
implementation  the  Poisson  equation  is  solved  with  the  Dirichlet  BCs,  and  the  resulting 
solution  is  constrained  by  the  Neumann  BCs.  The  Neumann  conditions  (a.k.a.  solvability 
constraints)  are  inserted  into  the  Helmholtz  problem  as  tau  conditions.  Starting  with  (2.8) 
where  the  Dirichlet  BCs  are  implemented  as  tau  conditions,  the  equation  has  the  form 

(2.12) 
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where  +  RP  -  J  is  the  postconditioned  Laplacian  with  the  t  conditions 

included.  Now  the  resulting  stream  function,  ,  is  forced  to  meet  the  Neumann 

conditions.  Since  the  Laplacian  with  Dirichlet  BCs  is  well  posed  and  invertible,  applying 
the  inverse  Laplacian  and  the  postconditioner  to  both  sides  of  (2.12)  yields 

(2.13)  f,  =  PV;  =  P"  (vj  )  ‘  ) 

the  Neumann  equations  are  implemented  by  defining  the  function 

fp  =  ^0,l,4,9,...,i^,...,M^^ .  Note  multiplication  of  a  Chebyshev  coefficient  vector  by 

is  identical  to  function  evaluation  of  the  derivative  at  the  boundary  r  =  1 ,  i.e. 

fp  •//„  =  .  Multiplying  into  (2.13)  and  applying  the  constraint  (2.10)  yields 

(2.14)  f.  ii?,  =  n  ■  P"  (V?  )■'  (-R^  ) = 0 

Define  the  row  vector  * ,  and  -R^  becomes  the  solvability  constraint 

on  vorticity  seen  by  applying  the  definition  to  (2.14) 

(2.15) 

To  find  the  solvability  constraint,  ,  multiply  its  definition  by  the  Laplacian  and 
transpose  the  resulting  equation  for 

(2.16)  (V^  ^  5/  =  (f,  •  f  forn  >  0 

Once  the  ’s  are  found  (2.15)  gives  the  exact  form  for  the  conditions  to  be  included  with 

the  Helmholtz  problem.  These  solvability  constraints  are  written  as  t  conditions  in 
equation  (2.9).  For  the  zero  Fourier  mode,  n  =  0 ,  the  Poisson  equation  simplifies  giving 


Integrating  over  the  disk  gives 

The  left  hand  side  simplifies  via  the  fundamental  theorem  of  calculus  and  matches 
condition  (2.10)  for  n  =  0 
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dVo 

dr 


r=l 


=  2nU 


The  right  hand  side  of  (2.17)  matches  the  definition  for  the  circulation 


(2.19) 


rlTT  plTl 

-Jo  ^^(OQvdrdO  =  j^cordrdO  =  -C 

Combining  the  results  of  (2.18)  and  (2.19)  equates  circulation  to  a  constant  value, 

C  =  -ItiU  .  Furthermore,  the  constant  circulation  implies  ^  =  0 .  Recall  from  section 
1.2.6  the  time  evolution  of  circulation  is  given  by 
(2.20)  f  =  +if^f^a>’rdr-fc 

and  the  resulting  value  for  the  circulation  is 


(2.21) 


r*  * 

C  =  e  ''  +  2^  ft)  rdr 
Jo 


Rendering  (2.21)  while  including  the  results  in  (2.19)  and  (2.18)  shows  the  tangential 
velocity  at  the  wall  is  a  function  of  the  initial  forcing 

C  =  2;r  f '  ft)*r<ir  =  -IttU 
Jo 

Incorporating  both  the  representations  for  circulation  into  (2.20)  gives 


=  -iTivd-COr, I  ,  +  Cco*rdr-^C  =  -2nv dMr. I  ,  - J 

dt  r  U|^=i  Jq  r  y)\r=\ 


-f/  =0 


This  produces  the  vorticity  constraint  for  Fourier  mode  zero.  The  condition  has  been 
transformed  from  the  stream  function  and  is  applied  to  the  vorticity  as  a  t  condition  in 
the  Helmholtz  equation 

9r^oU=0 

2.5.4  Accurate  initiation  of  the  time  integration  scheme 

The  3'^  order  SBDF  has  a  nominal  accuracy  of  the  order  j  where  is 

the  size  of  the  time  step.  In  order  to  keep  the  high  accuracy  during  the  initial  time  steps, 
the  solutions  are  initially  computed  with  a  order  SBDF  scheme  using  a  time  step  of 

{St'f .  Define  the  order  time  step  as 


6,={6t)^ 
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The  1®‘  order  scheme  needs  to  be  repeated  for  y,  steps  in  order  to  find  the  initial 

values  for  the  2"^  order  scheme.  To  start  the  2"^  order  SBDF  requires  the  initial 

2  2 

conditions  and  the  flow  parameters  at  time  {dt^  because  S2  ~  {St)  .  Note:  the  time  step 


here  is  only  approximate  since  ^  is  not  necessarily  an  integer.  The  actual  number  of  S^ 
time  steps  is  the  natural  number  j  such  that  -y<  +  Asa  result  the  time  step  for 

the2"%rderSBDFis  ' 

S2  =  jiStf 

To  find  the  initial  values  for  the  3*^^  order  SBDF,  the  2"^  order  scheme  runs  for  2  j 
iterations.  The  values  required  to  start  the  order  scheme  are  the  initial  value,  the  value 
at  time  ,  and  the  value  at  time  2 7^2  •  Thus  the  final  time  step  is 

(^3=/  {Stf^St 

The  error  between  the  actual  time  step  and  the  desired  time  step  will  equal  zero  if  St  is 
carefully  selected  so  that  ^  is  a  natural  number.  In  cases  where  is  not  a  natural 

number ,  the  maximum  time  step  error  with  7  <  -Jj-  + 1  is  given  as 

S2-St<{l  +  Stf  St-St  =  2{Stf  +{Stf 

Using  the  j  accurate  start  up  routine  adds  a  total  of  3j-2  steps  when  compared 

to  the  O  {St)  start  up. 


2.5.5  Time  integration  scheme 

An  implicit-explicit  time  integration  method  is  implemented  since  the  equations 
contain  both  linear  and  nonlinear  terms.  Linear  terms  and  time  integration  are  included 
implicitly  using  a  3'^^  order  BDF  scheme  while  the  nonlinear  terms  enter  explicitly, 
Ascher  etal[l].  The  combination  of  implicit  and  explicit  terms  provides  an  accurate 
calculation  with  a  large  stability  region.  The  3^^*  order  SBDF  has  the  form 


(2.22) 


-3or* ^ ^ 

^ - ^ - 2 - =  3/?^  -  3/?*"‘  + 


St 


where  a  is  calculated  implicitly,  P  contains  the  nonlinear  terms  and  is  calculated 
explicitly,  St  is  the  time  step,  and  k  is  the  time  level.  For  our  problem  the  Poisson 
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equation  is  solved  explicitly  at  the  current  time  level.  The  vorticity  at  a  given  time 
produces  the  stream  function  at  the  same  time  level 

(2.23)  VV*  = 

The  inverse  Laplacian  as  applied  to  the  vorticity  is  resolved  using  the  postconditioning 
method  discussed  later.  Once  the  stream  function  and  vorticity  are  on  the  same  time 
level,  the  Jacobian  is  calculated  and  the  Helmholtz  equation  produces  the  next  vorticity 
via  the  SBDF  method.  The  exact  form  for  the  calculation  is 

(2.24) 

-x)j£y^+i  +<^^(3/  + 

where  /*  =  ^co*  -  J {0/  ,\i/^ ) .  At  this  point  the  vorticity  is  recast  with  the  multiplier 

in  order  to  smooth  the  flow  near  the  origin  singularity.  In  general  the  vorticity  is 
represented  as 

(o{r)  =  r^a)°{r) 

and  the  Helmholtz  equation  is  solved  for  co°(r) .  Variations  of  this  representation  exist 
when  solving  Fourier  modes  0  or  1,  and  a  full  description  is  provided  later. 

2.5.6  Diagnostics 

The  primary  diagnostics  for  the  code  are  completed  using  the  energy  evolution, 
enstrophy  evolution,  and  angular  momentum  evolution.  Recall  from  section  1  the  energy 
evolution  is  given  by 

(2.25)  ,d0-vW+^-n^''u-u*d0rdr-i^T 
the  enstrophy  evolution  by 

^  («'  X  -  2vJ'7o  +  {}de^fydrd0  + 

(2.26) 

and  the  angular  momentum  evolution  by 

(2.27)  ^  =  -2OT'ffl„  (1) + 4w»„  (1)  -  -^1 V"  [»;  (r)  -  {r)]dr 

To  check  the  accuracy  of  the  code  the  left  hand  sides  of  (2.25),  (2.26),  and  (2.27)  are 
computed  over  three  time  steps  using  a  second  order  difference  scheme.  The  numerically 
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computed  derivatives  are  subtracted  from  their  respective  right  hand  sides,  and  the 
relative  error  is  computed. 

2.6  Forced  Periodicity 

Investigating  secondary  bifurcations  from  steady,  inhomogeneous  flow  requires 
the  classification  of  stable  solutions  for  various  values  of  the  Reynolds  number.  As  with 
the  initial  bifurcation  from  the  steady  homogeneous  solution,  the  Reynolds  number  is  the 
appropriate  bifurcation  parameter.  To  find  the  bifurcation,  the  steady  solution  is  forced 
into  an  unstable  regime.  Artificial  stability  of  a  flow  is  realized  by  forcing  the  periodicity 
of  the  stable  solution  beyond  the  critical  value  of  the  Reynolds  number.  The  periodicity 
of  the  solution  is  enforced  by  equating  all  other  Fourier  modes  to  zero  e.g.  a  period  4 
solution  is  forced  by  allowing  only  Ak  Fourier  modes  where  0<k< .  Using  the 

continuous  evolution  of  eigenvalues  through  a  change  in  Reynolds  number,  stable 
solutions  can  be  tracked  onto  their  associated  instability  regions. 

2.7  Conditioning  of  the  problem 

The  conditioning  for  the  numerical  simulation  follows  directly  from  the 
Helmholtz  conditioning,  Torres  &  Coutsias  [32].  The  Poisson  equation  is  well  posed  and 
robust  to  small  perturbation.  The  Helmholtz  equation  is  derived  by  applying  the  SBDF 
scheme  to  equation  (2.9).  The  resulting  equation  is 

(2.28)  +  5RI[,]P  +  (4  -  n"  )p"  +  il]®;  =  I, 

where  has  the  same  form  as  the  right  hand  side  of  (2.24).  In  general  the  linear 
Helmholtz  equation  has  the  form 

(2.29)  Aft)°  =  g 

where  cd°  is  used  to  find  the  vorticity  through  the  smoothing  and  postconditioning 

(2.30)  €0  =  for  ;  =  1, 2 

To  estimate  the  condition  of  the  problem  consider  a  perturbed  linear  system  (2.29)  and  a 
perturbed  conversion  (2.30) 

(2.31)  (A  +  <5A)(d)°+(5d)°)=g  +  ^g 

(2.32)  Q)  +  d(d  =  (r^  +  ^R^'  )(P  +  ^P)(d)“  +  Sa>° ) 
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The  conditioning  for  the  problem  is  measured  through  the  magnitude  of  the  relative  error 
in  vorticity.  This  error  compares  the  norm  of  the  vorticity  error  to  the  norm  of  the 
vorticity 

IS  =  condition  number 

HI 

The  details  for  the  formulation  of  the  condition  number  are  given  in  appendix  B,  but  the 
result  is  stated  here  for  convenience 

(2.33)  M  <  M|R;p2  A-'  I  (;|r(g)  +  22r(  A)) 

where  A  +  5W(,jP  +  (4-n^)p^ +^l],  0<a<tt,and 
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3  Stability  Analysis 

3.1  Introduction 

The  fundamental  concept  for  stability  analysis  addresses  how  a  physical  system 
evolves  in  time  when  considering  conditions  near  one  of  the  system’s  equilibrium  states. 
The  flow  is  perturbed  slightly  about  the  state  in  question,  and  the  evolution  of  the 
perturbation  is  studied.  When  only  terms  linear  in  the  perturbation  are  considered 
(linearized  analysis),  results  about  the  stability  of  the  given  state  are  obtainable,  but  in  the 
case  of  instability  no  conclusion  can  be  drawn  about  the  eventual  state  that  will  be 
reached.  From  experimental  work,  Rabaud  &  Couder  [25],  Niino  &  Misawa  [23],  and 
Fruh  &  Read  [15],  the  bifurcation  from  the  homogeneous  flow  is  known  to  be  a 
supercritical  pitchfork.  To  derive  the  form  of  the  bifurcation  analytically  one  must 
extend  the  treatment  to  a  weakly  nonlinear  analysis  which  includes  terms  that  can 
produce  saturation  of  the  growing  perturbation.  Such  analysis  for  2D,  rotating  fluids  has 
been  carried  out  by  Bergeron  et  al  [3],  Churilov  &  Shukhman  [7],  and  Niino  &  Misawa 
[23]  confirming  the  experimental  findings  for  the  first  transition  from  rotational 
symmetry  to  small  amplitude  rotating  waves.  In  this  stage  of  the  work  we  limit  our 
analysis  to  the  linearized  problem. 

Initially,  a  specific  description  of  linear  stability  is  provided.  Beginning  with  an 
equilibrium  state,  a  perturbation  is  introduced  so  that  the  robustness  of  the  equilibrium 
solution  can  be  investigated.  With  e  as  the  perturbation  parameter,  the  perturbed 
solution  resides  within  an  e  -  ball  of  the  equilibrium  solution,  see  figure  3.1.  The  steady 
system  under  investigation  is  autonomous  implying  that  the  temporal  dependence  for 

disturbances  from  the  equilibrium  is  of  the  form  ,  where  AeC .  Evolution  times 
associated  with  the  unstable  solutions  are  finite  since  perturbed  solutions  are  only 
eonsidered  within  the  £  -  ball .  This  technique  addresses  the  initial  stability  of  the 
solution  for  a  given  perturbation. 
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fig  3.1  Linearized  solution  and  bifurcation. 

The  limitations  of  this  method  need  to  be  addressed.  As  previously  stated, 
unstable  solutions  cannot  be  considered  for  times  beyond  which  the  solution  is  no  longer 
contained  within  the  s  -  ball .  Saturation  of  the  unstable  solution  may  occur  beyond  the 
limits  of  this  analysis.  However,  the  nonlinear  terms  may  drive  the  perturbed  solution 
away  from  a  solution  found  to  be  stable  via  the  linear  stability  analysis.  This  second 
issue  does  not  pose  a  problem  in  our  case  since  the  spectrum  for  this  system  is  bounded 
away  from  the  imaginary  axis,  and  for  such  operators  the  linear  analysis  is  sufficient, 
Trefethen  [xx].  In  addition,  direct  numerical  simulation  demonstrates  that  the  linear 
behavior  dominates  the  flow  evolution  near  the  bifurcation. 

3.2  Homogeneous  Flow 

The  homogeneous  flow  is  characterized  via  the  fact  that  the  velocity  has  no 
variations  in  the  6  direction.  The  flow  is  steady  in  time,  and  the  perturbation  is 
introduced  into  both  the  stream  function  and  vorticity.  The  problem  can  be  uncoupled  in 
the  Fourier  expansion  so  as  to  consider  only  a  single  Fourier  mode  when  applying  a 
perturbation  to  the  problem.  To  analyze  this  problem  we  consider  the  dimensionless 
vorticity-stream  function  equations 

(3.1)  a,o  +  7(fy,V^)  =  -^[v^(y  +  -^(ty*-ft))] 

(3.2)  vV=-fy 

where  Re  is  the  Reynolds  number.  Re  =  is  the  angular  velocity  of  the  rigidly 

rotating  flow  inside  radius  a;  is  the  characteristic  length  scale;  and  (O*  is  the  force 
applied  as  a  vorticity  field. 

For  the  linear  stability  analysis  solutions  are  perturbed  from  steady  state.  Our 
interest  lies  in  time  evolution  of  the  instabilities  in  an  autonomous  system.  Consider  the 


steady  state  problem  and  write  the  Laplacian  operator  for  polar  coordinates  as 
'^^^0  ■  The  resulting  system  of  equations  from  (3.1)  and  (3.2)  is 

(3.3)  -X^rO^eW  -  ^rV^eCo)  =  +  +  ^ (of  -  ty)] 

(3.4)  dly/+^d,y/+jr^ly/  =  -co 

Using  the  system  (3.3)/(3.4),  we  consider  the  stability  of  the  azimuthally 
symmetric  flow.  For  the  equilibrium  state  (O  and  y/  are  only  functions  of  r ,  and  we 

write  the  vorticity  and  stream  function  as  tyo(r)  and  y/^ir) ,  respectively.  Due  to  the 
simplicity  of  the  base  state,  we  consider  perturbations  of  the  form  ... 

(3.5)  ty  =  tyo(r)  +  f<w„;i(r)e^'^"'^ 

(3.6)  y/  =  y/(^{r)  +  sy/,,.^{r)Q^'^^’^^ 

where  e  is  an  amplitude  parameter,  n  is  the  number  of  the  Fourier  mode,  and  A  is  an 
eigenvalue  of  the  n*  mode.  We  seek  the  form  of  the  corrections  co^-iir)  and  y/„.iir)  as 
well  as  the  corresponding  eigenvalues. 

The  6  independence  of  the  base  state  allows  uncoupling  of  the  Fourier  modes 
and  provides  for  the  independent  consideration  of  each  Fourier  mode.  After  allowing  the 
substitution  of  y  =  At-\-  inO ,  the  asymptotic  expansions,  (3.5)  and  (3.6),  are  substituted 
into  equations  (3.3)  and  (3.4).  The  time  dependency  of  the  perturbation  forces  the 
inclusion  of  the  time  derivative  term,  .  The  resulting  system  is: 

(3.7) 

e(0„.,M  +\\d,o)Q  (einyr„.^t^)-d,yfQ  )]  =  )+ 

(3.8)  dVo  +  +7(^r^«^o  +  " “^o  “ 

Analysis  beyond  this  point  considers  the  system  as  two  separate  pieces  at  distinct  orders. 
The  order  one,  O  (l),  terms  are  grouped  separately  from  the  order  epsilon,  O  (e),  terms. 

The  O  (l)  system  involves  only  the  homogeneous  flow  without  perturbations.  For 
Reynolds  numbers  below  the  critical  value,  the  O  (l)  equations 


48 


(3.9) 


(3.10)  d]\i/Q+\d,iifQ=-a)Q 

give  the  same  stream  function  and  vorticity  information  as  the  numerical  simulation. 

The  )  equations  pertain  to  the  linear  perturbation  of  the  homogeneous 

problem.  Since  the  system  is  linear,  eigenvalue  decomposition  can  be  used  to  give  the 
stability  of  the  flow  for  various  Reynolds  numbers.  A  positive  eigenvalue  corresponds  to 
an  instability  in  the  homogeneous  flow,  and  during  our  investigation  we  search  for  a  pair 
of  eigenvalues  that  cross  the  imaginary  axis.  The  O  (e)  systems 


(3.11) 

+  y  A;0  - 

(3.12) 

are  coupled  through  the  vorticity  and  stream  function  perturbations.  Using  the 
homogeneous  flow  parameters  found  by  solving  (3.9)  and  (3.10),  the  system  (3.11)/(3.12) 
is  reconstituted  as  a  generalized  eigenvalue  problem,  .  For  equation  (3.11) 

the  partial  derivative,  and  .  ^re  required  once  the  homogeneous  vorticity 

and  stream  function  are  known.  In  order  to  find  these  terms  directly,  (3.9)  and  (3.10)  are 
differentiated  with  respect  to  r .  Applying  to  the  homogeneous  equations  yields 

(3.13)  (aj  +13, 

(3.14)  (aj  +13,  -y)a,(i?„i, = -a,a..„ 


where  the  forcing  is  identically  zero  for  all  nonzero  Fourier  modes,  d)*  =  0  for  n  >  0 . 
Thus  the  stability  problem  is  resolved  by  solving  (3.13)  and  (3.14)  for  d;.d)„.o  and  , 
and  formulating  the  generalized  eigenvalue  problem  from  (3.11)  and  (3.12). 

To  begin  the  singularities  are  removed  through  a  multiplication  by  .  Recall  the 
functions  ,  and  are  functions  of  r  only,  not  0 ,  implying  they  have 

only  Chebyshev  expansions.  Each  of  the  functions  above  has  an  expansion  via 
Chebyshev  polynomials  e.g.  d)„_o  =  [^y„,o;o  ^n.i,o  •" 
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Wn,\=\Wnsy>\  Vn,\\\  Thc  entry  in  thc  vcctor  is  thc  Coefficient  for 

the  Chebyshev  polynomial,  namely  Tyfr) .  With  the  Chebyshev  vector 

representation  of  vorticity  and  stream  function,  multiplication  by  r  and  differentiation  by 
r  are  accomplished  using  matrix-vector  multiplication  as  described  in  appendix  D. 

Using  the  Chebyshev  coefficient  vectors  in  combination  with  R  and  D  matrices  the 
equations  (3.13)  and  (3.14)  are  written  in  matrix-vector  form.  After  applying  the 

multiplication  by  r  the  0(1)  equations  are 

(3.15)  (rV  +RD-I-iR^)3A,o 

(3.16)  (r"D"  +RD-l)a,(i?„.„  =-R"3A.o 


These  equations  are  straight  forward  to  solve  since  d^ca*  is  a  given  function 

3A=^sech^(^)[^|tta„h(^y 

The  resulting  solutions  to  (3.15)  and  (3.16)  are  the  vectors  dr^„  o  and  d^^„  Q ,  and  system 
(3.1 1)/(3.12)  can  now  be  written  in  matrix-vector  form  and  solved. 

To  this  end  apply  the  multiplication  by  and  change  system  (3.1 1)/(3.12)  into 
matrix-vector  representation 


(dr^„;0  f  -  drl^„,c 

(3.17) 

+  RD  -  n^I -J-R2  )ft), 

(3.18) 

(r^D^  +  RD  -  =  -R^d)„., 

where  *  represents  the  convolution  of  two  vectors  in  Chebyshev  mode  space.  The 

convolution  terms,  and  *  ^n;i  >  rewritten  as  a  matrix-vector  product 

where  and  are  the  unknowns  and  retain  their  vector  form.  As  a  result  the 

and  vectors  are  written  as  Chebyshev  convolution  matrices.  The  details  for  writing 

a  vector  as  a  Chebyshev  convolution  matrix  are  given  in  appendix  E.  Rewriting  (3.17) 
and  (3.18)  while  including  the  Chebyshev  convolution  matrices  yields 
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Now  to  find  a  specific  solution  to  the  generalized  eigenvalue  system,  boundary 
conditions  (BCs)  are  required.  The  physical  BCs  for  this  flow  problem  are  already 
satisfied  through  the  O  (l)  equations.  The  BCs  for  the  O  (e)  equations  are  homogeneous 


since  the  perturbation  has  homogeneous  end  conditions.  For  the  no-slip  flow  at  the 
boundary,  Dirichlet  and  Neumann  conditions  are  applied  to  the  stream  function 


(3.22) 


^rl^„U±l=0  for  n>0 

3.roU.=£/ 


(3.23)  <>'.^,=0 


The  Neumann  conditions  are  transferred  to  the  vorticity  as  solvability  constraints  as 
explained  in  section  2.5.3.  The  Dirichlet  conditions  are  enforced  at  both  ends  of  the 
radial  grid.  We  note  that  all  Chebyshev  polynomials  evaluate  to  1  at  r  =  1 , 
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y/  >  0;  r^l  _j  =  1 ,  and  the  Chebyshev  polynomials  are  equal  to  ±1  at  r  =  -1 , 

V;  >  0;T'y|  _  =  (-ly .  The  vector  representations  for  the  Dirichlet  BCs  are 

0  =  [1  1  •••  lf@r  =  l 
0  =  [1  -1  1  -1  •••  if  @r  =  -l 

These  conditions  replace  the  last  two  stream  function  equations  on  both  sides  of  (3.21). 
The  solvability  constraints  for  the  vorticity  are  the  vectors,  ,  found  by  solving 

These  constraints  replace  the  final  vorticity  equations  on  both  sides  of  (3.21). 

3.3  Solving  the  generalized  eigenvalue  problem, 

Several  iterative  methods  are  available  to  solve  this  problem.  Due  to  the 
relatively  small  size  of  these  computations,  the  generalized  eigenvalue  problem  is  solved 
using  MATLAB’s  QZ  implementation.  The  eigenvalue  problem  requires  the  O  (l) 

vorticity  and  stream  function,  and  equations  (3.15)  and  (3.16)  are  solved  first  for  the 
vectors  ^nd  drW„.fi  ■  With  these  initial  flow  parameters,  system  (3.21)  is  solved 

and  the  spectrum  is  recaptured  by  removing  the  Ekman  coefficient  term 
A  =  ju - 

Each  Fourier  mode  is  solved  separately  so  the  issue  of  resolution  relates  only  to  the 
Chebyshev  modes.  The  generalized  spectrum  contains  the  dynamics  of  the  physical  flow 
as  well  as  spurious  eigenvalues  associated  with  the  numerical  implementation  of  the 
operators.  Consider  the  problem  solved  with  a  variety  of  Chebyshev  modes,  figure  3.2. 

In  order  to  separate  the  physical  spectrum  from  the  pseudo-spectrum,  a  perturbation  is 
introduced  to  the  problem.  To  perturb  the  generalized  eigenvalue  problem,  complete 
spectra  are  found  for  three  distinct  Chebyshev  resolutions,  and  the  eigenvalues  are 
compared.  The  eigenvalues  that  differ  only  by  the  magnitude  of  the  error  resulting  from 
the  different  truncations  are  kept  as  these  represent  the  physical  spectrum.  The  “moving” 
eigenvalues  are  part  of  the  pseudo-spectra  and  are  discarded. 
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fig  3.2  Eigenvalue  plots  demonstrating  increasing  resolution  of  the  physical  spectrum 
corresponding  to  an  increasing  number  of  Chebyshev  modes. 

The  number  of  computations  required  to  solve  the  problem  is  quartered  by 

noticing  that  equation  (3.21)  is  even.  The  even  nature  of  the  problem  implies  only  even 

products  of  Chebyshev  and  Fourier  modes  exist.  The  non  zero  elements  of  the 

(M  +1)x(M  + 1)  Chebyshev  matrices  appear  in  the  following  checkerboard  pattern 
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Due  to  this  structure  the  boundary  conditions  are  rewritten  as  linear  combinations  of  one 
another  thus  matching  the  even  or  odd  parity  of  the  problem. 

After  the  BCs  have  been  placed  in  the  generalized  eigenvalue  problem,  the 
stability  of  each  Fourier  mode  can  be  determined  by  the  resulting  eigenvalues.  An 
unstable  Fourier  mode  indicates  the  flow  has  bifurcated  to  a  new  solution. 

3.4  Secondary  bifurcations 

Here  we  are  looking  for  secondary  bifurcations  from  a  steady  state,  inhomogeneous 
solution.  The  steady  solution  differs  from  the  homogeneous  analysis  since  the  solution 
varies  in  6  as  well  as  r ,  and  is  characterized  by  a  specific  number  of  vortices.  Since  the 
flow  is  now  dependent  upon  azimuthal  variations,  the  analytical  perturbations  from 
steady  flow  are  functions  of  r ,  6  ,  and  t.  As  in  the  homogeneous  problem,  only  terms 
linear  in  the  perturbation  are  considered.  The  forms  for  the  vorticity  and  stream  function 
steady  solutions  with  perturbations  are 

(3.24)  co(r,0,t)  =  cOQ(r,0)  +  sa)i(r,0)e^'^ 

(3.25)  }//ir,0,t)  =  i//Q(r,0)  +  si/r^ir,0)e^' 

Unlike  the  homogeneous  stability  analysis,  the  steady  flow  characterized  through  cOq  and 
y/f)  is  a  function  of  0  and  is  obtained  from  the  numerical  simulation.  Of  interest  for  this 
perturbation  study  are  the  0(e)  amplitudes  (O^  and  y/i,  and  their  time  evolution  given  by 

A .  In  the  generalized  eigenvalue  framework  A  is  the  eigenvalue  and  ty,  and  y/^  are  the 
eigenfunctions  for  this  problem.  Unlike  the  homogeneous  problem,  the  system  of 
equations  do  not  uncouple  through  the  Fourier  modes  and  must  be  solved  as  a  system. 

Our  interest  lies  in  time  evolution  of  the  instabilities  in  an  autonomous  system. 
The  resulting  system  of  equations  from  (3.1)  and  (3.2)  is 

(3.26)  ^{d.codgy/  -  d.yrdgO))  =  +  ty)] 

(3.27)  dly^  +  ^d,y/  +  jrdly/  =  -a) 

Applying  the  perturbed  solutions  from  (3.24)  and  (3.25)  to  the  system  (3.26)7(3.27)  gives 
the  following  system  of  equations: 
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(3.28) 


eAcOiC^'  +j(d,o)od0i/ro  -d.i/^od^coo  +  ed.coidffi^o^^^)- 
^(ed.i/rodffCOie^'  +  )  = 

■^[v^cyo  +  +^((o  -(o^- eco^q^  )] 

(3 .29)  =  -^0  ~ 

This  perturbed  system  is  time  independent  and  can  be  separated  according  to  two  distinct 
magnitudes,  0(1)  and  0(e).  Here  (Oq  and  y/Q  are  already  known  via  numerical 

simulation  of  the  flow  so  the  0(l)  equations  are  used  only  to  find  d^cb^  and  d^y/Q 
directly  instead  of  resorting  to  matrix  vector  multiplication,  namely  Dft)o  and  Y)y/Q . 

From  the  system  (3.28)/(3.29)  the  0(l)  equations  are 

(3.30)  -«o)] 

(3.31)  vVo  =  -y>o 

The  primary  equations  considered  for  stability  analysis  are  the  0{e)  equations 

+):{d,(0^dgyfQ&^'  -d^y/^dgCo^e.^'  +d^o)f^dgy/^c^'  -  d^y/^dgCOff.^‘')  = 

(3.32)  '  , 

(3.33)  vVie^'=-y>ie'^' 

which  contain  a  variety  of  nonlinear  terms  resulting  from  the  Jacobian.  Using  the  matrix- 
vector  technique  again,  system  (3.32)7(3.33)  is  converted  to  vector  form  through  the 
Chebyshev  and  Fourier  modal  coefficients.  Here  the  form  differs  from  the  homogeneous 
problem  since  coupled  Fourier  modes  are  included  in  the  steady  flow,  and  the  coupled 
modes  negate  the  approached  of  separating  the  Fourier  modes.  A  new  matrix- vector 
form  is  required  where  the  Chebyshev  and  Fourier  modes  are  included  in  a  single  vector. 
The  following  expansion  is  generic  and  is  used  for  both  the  stream  function  and  vorticity 

2  M 

(3.34)  C (^^)  =  E  E  (Cm  cosn0  +  Q,„  sin nd)T^  (r) 

n=0m=0 

where  ^  is  either  stream  function  or  vorticity  ( or  (O),  T,,,  (r)  is  the  m'*  Chebyshev 
polynomial,  and  ^nm  the  ^  cosine  and  sine  coefficients,  respectively,  for  the  n'* 
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Fourier  mode  and  the  m'*  Chebyshev  mode,  and  the  domain  for  (r,d)  is  -1  <  r  <  1 ; 

0<d<27r .  This  expansion  establishes  the  Chebyshev  coefficients  as  the  inner 
expansion  i.e.  for  each  Fourier  mode  the  Chebyshev  coefficients  are  together.  An 
alternate  characterization  is  used  when  only  considering  the  outer,  Fourier,  expansion: 

(3.35)  C (a^)=  E(C  {r)cosn9  +Q  (r)sinn(9) 

n=0 

where  ^  is  either  stream  function  or  vorticity  (y/  or  co)  and  (r)  and  Cn  (^)  the 

Chebyshev  expansion  for  the  Fourier  cosine  and  sine  modes.  Translating  the  above 
expression  into  a  vector  of  coefficients  shows  the  Chebyshev  expansion  as  the  inner 
expansion  and  the  Fourier  expansion  as  the  outer  expansion.  Considering  (3.35)  in  full 
vector  form  produces  the  Fourier  vector 

f(r,e)  =  [f„'(r)  f,'(r)  f|(r)  0  ff(r)  ft(r)  -  f|(r)f 

Thus  each  element  in  the  Fourier  vector  is  a  Chebyshev  coefficient  vector  having  the 
same  form  as  the  one  described  for  the  homogeneous  problem,  and  the  result  is  a  single 
column  vector  containing  {M  +  l){N  +  2)  elements.  Note  the  sine  zero  coefficients  (all 

zeros)  are  kept  to  maintain  symmetry  in  the  Fourier  vector. 

With  the  matrix-vector  form  defined  we  consider  solutions  to  the  system  (3.32) 
/(3.33).  Once  again  elimination  of  the  singularity  occurring  at  r  =  0  is  the  first  point  of 

consideration.  Multiplication  by  will  remove  the  point  singularity.  To  maneuver 
within  the  full  vector  space  defined  above,  operations  such  as  differentiation  by  r  and 
differentiation  by  9  must  be  rendered  in  the  appropriate  matrices.  As  with  the 
homogeneous  problem,  operations  on  the  Chebyshev  basis  are  independent  of  9 ,  and  this 
implies  the  matrices  for  multiplication  and  differentiation  by  r  are  block  diagonal 

R  =  I,,3®R 

D  =  I,,,0D 

where  is  the  AT  -1-2  xiV  +  2  identity  matrix.  Within  the  trigonometric  basis 
differentiation  by  9  corresponds  to  reassignment  of  the  cosine  and  sine  coefficients 
while  scaling  by  the  mode  number 
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0 

where  increments  as  the  n  - 1  row  of  .  Detailed  derivations  of  these  matrices 

can  be  found  in  appendix  D. 

Now  that  all  the  tools  are  in  place,  rewrite  (3.32)  and  (3.33)  in  matrix-vector  form. 
First  consider  the  Laplacian  operator  for  polar  coordinates,  +79,.  +-;j^0  •  The 

point  singularity  is  removed  via  multiplication  by  ,  and  after  applying  the  previous 

Ak  ys  ^  /\  ys  Ak  ^ 

definitions,  r  V  =  R  D  -I-  RD  -i-  Q  .  Using  this  representation  of  the  Laplacian  and 
converting  the  remaining  portions  of  the  equations,  gives  the  system  in  matrix-vector 

AR^Q)^  -l-RQ|^o*Dd),  -l^^^o*Qd),  -l-RDd)o*Q^]  -RQd)o*D^«^j  = 
-^(R^D^d)i  -f-  RDd)i  +  Q2ft)i  -  R^  ) 

R^D^^j  +  to  =  -R^d)i 

where  *  is  convolution  in  Fourier  space.  This  system  is  similar  to  the  homogeneous 

F 

yv 

problem  except  for  the  addition  of  the  Q  matrix.  The  existence  of  coupled  Fourier 
modes  increases  the  number  of  nonlinear  terms  resulting  from  the  Jacobian.  These 
nonlinear  terms  in  point  space  translate  to  convolutions  in  Fourier  mode  space, 

RQ^q*D<Wj,  — RD^Q*Q^yj,  RDc^q^Q^j  ,  and  -RQd>o*D^^i.  The  convolutions  are 

resolved  by  transforming  the  leading  vector  into  its  associated  convolution  matrix.  Since 
the  generalized  eigenvalue  calculation  is  not  iterative  with  respect  to  the  Fourier 
convolutions,  the  use  of  fast  transforms  such  as  the  FFT  is  not  required. 

The  Fourier  convolution  matrices,  RQ^i^o®,  RDd)oQ,and  -RQtyoD,are 

formed  following  the  detailed  procedure  in  appendix  F. 

The  formulation  of  the  generalized  eigenvalue  problem,  =  juC^ ,  can  be 
accomplished  in  several  ways.  With  both  the  vorticity  and  stream  function  as  variables, 
system  (3.36)/(3.37)  is  written  as 


form 

(3.36) 

(3.37) 


57 


(3.38) 


+  RD  +  Q^ )-RQt^oD  +  RD^oQ  RQ^o® " RD^oQ 


R' 


R^D^+RD  +  Q" 


R^  0 
0  0 


where  ^  = 


0)^ 


and  ^  =  X->r  — ^- .  The  size  of  the  matrices  in  (3.38)  is 

Re/i 


2(M  + 1)(//  +  2)x2(M  +  1)(N  +  2)  and  this  size  of  matrix  resolves  a  problem 


containing  0-^M  Chebyshev  modes  and  0  — >  y  Fourier  modes.  To  reduce  the  size  of 


the  impending  eigenvalue  calculation,  system  (336)1(331)  can  be  written  in  terms  of 
only  the  vorticity.  Instabilities  in  the  flow  are  captured  as  accurately  when  considering 
only  the  vorticity  since  the  onset  of  the  instability  is  the  only  feature  of  interest.  This 
analysis  does  not  make  any  claims  about  the  flow  after  an  instability  has  occurred.  To 
write  the  problem  in  terms  of  vorticity  only  begin  by  rewriting  (3.37)  as 

R^D^+RD  +  Q^J  R^d)i 


The  inverse  Laplacian  operator  is  defined  when  boundary  conditions  are  applied,  and  this 
provides  a  well  posed  problem.  For  the  no-slip  flow  Dirichlet  boundary  conditions  are 
applied  to  the  stream  function.  Once  the  stream  function  is  written  in  terms  of  the 
vorticity,  equation  (3.39)  is  substituted  into  (3.36)  resulting  in  a  single  vorticity  equation 

/S  .A.  .A  /  A  A  A  AA  AA  ^  O 

RQt^oD«)i  +  (RQft>oD  -  RD^oQ)[R  +  RD  +  ]  R^w,  - 


(3.40) 


-  (R^D^  +  RD  +  jd),  -  RD 


In  the  pure  vorticity  form  only  the  convolution  terms  contain  non  trivial  blocks  off  the 
main  diagonal. 

Boundary  conditions  (BCs)  are  the  final  requirement  for  solving  this  problem.  As 
before  with  the  homogeneous  stability  analysis,  the  BCs  for  the  0(e)  equations  are 


homogeneous.  The  stream  function  Dirichlet  conditions  are  already  included  in  the 
Laplacian  of  (3.39).  These  conditions  are  applied  by  removing  the  smallest  magnitude 
equations  prior  to  the  inversion  of  the  Laplacian.  The  operator  is  block  diagonal,  and  the 
BCs  are  applied  by  replacing  the  last  two  rows  in  each  (M  +  l)x  (M  + 1)  block.  For  the 
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pure  vorticity  equation,  (3.40),  the  vorticity  solvability  constraints,  ,  are  enforced.  The 
solvability  constraints  are  determined  as  before 

and  enter  as  the  final  two  rows  in  each  (M  +  l)x  (M  + 1)  block  along  the  diagonal.  In 
addition  elements  in  the  off  diagonal  columns  are  required  to  be  zero. 

3.5  Solving  the  generalized  eigenvalue  problem, 

The  generalized  eigenvalue  problem  with  the  full  vorticity  vector  can  no  longer  be 
solved  with  the  MATLAB  QZ  solver.  The  size  of  the  problem  is  prohibitive  and  the 
LAPACK  F77  solver  (dggev.f)  is  used  instead.  The  block  structure  inherent  in  A  and  C 
allows  for  a  reduction  in  the  size  of  the  QZ  computation.  As  with  the  homogeneous  case, 
this  problem  is  even.  The  blocks  along  the  diagonal  are  even  while  the  blocks  off  the 
diagonal  match  the  parity  of  the  associated  Fourier  mode.  The  matrix  with  boundary 
conditions  has  the  form 
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The  form  of  the  vorticity  vector  also  alternates  parity  as  the  Fourier  mode 

d)i=[*  0  *  0  *  0  •••  0  0  0  0  0  •••  0  0  *  0  *  0  O]^ 
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The  vector,  A.O)^ ,  resulting  from  the  checkerboard  pattern  in  A  contains  only  — 

non  zero  values.  This  effectively  reduces  the  size  of  the  generalized  eigenvalue  problem 
by  a  factor  of  four.  The  only  relevant  terms  in  A  are 

A  [l :  2 :  (M  +  l)(A  +  2),1 : 2 :  (M  +  l)(iV  +  2)] 

and  the  actual  size  of  the  generalized  eigenvalue  problem  is  x 

Even  with  the  reduction  in  size,  the  primary  concern  for  this  problem  is  the  fact 
that  A  and  C  are  large.  Investigation  of  500  Chebyshev  modes  and  200  Fourier  modes 
requires  the  solution  of  a  2,500,000,000  element  generalized  eigenvalue  problem.  This 
severely  limits  the  size  of  the  Chebyshev-Fourier  expansion  that  can  be  considered. 
However,  two  facts  allow  this  generalized  eigenvalue  problem  to  be  solved  in  an  efficient 
manner.  Recall  the  {M  +1)x(M  +1)  Chebyshev  convolution  matrices  are  formed  from 

vectors  containing  2M  + 1  modes.  This  allows  for  the  consideration  of  500  Chebyshev 
modes  using  a  matrix  of  size  250  x  250.  Next,  the  forced  periodicity  of  the  solution 
allows  for  the  consideration  of  a  greater  number  of  Fourier  modes.  The  flow  field  in  the 
numerical  simulation  is  forced  into  a  specified  symmetry.  When  an  n  =  6  braid  is  forced 
only  the  n  =  6k  Fourier  modes  need  to  be  considered.  As  a  result  500  Chebyshev  modes 
and  200  Fourier  modes  form  a  problem  with  only  18,000,000  elements. 
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4  Results 

4.1  Introduction 

The  stability  of  a  two-dimensional,  rotating  fluid  in  a  bounded  geometry  has  been 
studied  through  experimentation  and  analysis  using  two  basic  geometries,  thin  annular 
regions  or  thin  cylinders  (disks).  The  disk  and  annular  geometries  exhibit  similar 
dynamics  for  low  Reynolds  numbers  as  shown  is  section  4.2.  The  initial  bifurcation  is 
robust  to  several  forms  of  forcing  and  varying  aspect  ratios'.  For  higher  Reynolds 
numbers  the  post  on  the  inside  of  the  annulus  interferes  with  the  evolution  of  the 
asymmetric  tripole  which  appears  to  be  the  high  Reynolds  number  limit  for  the  disk  flow 
as  shown  in  section  4.3.  The  Ekman  pumping  we  have  selected  closely  resembles  the 
forcing  implemented  in  the  experiments  by  Friih  &  Read  [15],  and  direct  comparisons  are 
made  between  their  observations  and  our  simulations.  For  transitions  beyond  the  initial 
bifurcation,  section  4.4  discusses  symmetry  breaking  bifurcations  leading  to  periodic  and 
irregular  states. 

An  adequate  parameterization  for  the  transition  between  states  is  provided  by  the 
Reynolds  number,  Ehde  &  Titman  [18].  In  general,  the  Reynolds  number  is  the  ratio  of 
inertial  force  to  viscous  force 

Re  =  -^ 

where  U  is  the  characteristic  velocity  scale,  D  is  the  characteristic  length  scale,  and  v  is 
the  viscosity.  Throughout  the  field  of  rotating  shear  layer  stability,  two  methods  are  used 
to  define  the  Reynolds  number,  Dolzhanskii  et  al  [11].  The  definitions  differ  in  the 
choice  of  the  characteristic  length  scale;  one  method  uses  the  cell  height  while  the  other 
uses  the  shear  layer  thickness.  Dolzhanskii  [10]  shows  the  quasi-two-dimensional 
problem  is  accurately  represented  by  using  the  shear  layer  thickness  as  the  characteristic 
length  scale.  Our  definition  for  the  Reynolds  number  is 

where  AQ  is  the  angular  velocity,  a  is  the  distance  to  the  shear  layer,  and  L  is  the  shear 
layer  thickness. 

*  aspect  ratio  is  defined  as  the  cell  height  divided  by  the  location  of  the  shear  layer 


4.2  Homogeneous  instability 

The  initial  onset  of  flow  instability  in  the  disk  mirrors  the  results  found  for  the 
annular  geometry,  Bergeron  et  al  [3].  Near  the  center  of  the  disk,  the  flow  is  similar  to 
rigid  body  rotation  which  produces  little  impact  on  the  initial  instability  when  compared 
to  the  annular  geometry.  To  demonstrate  this  idea,  we  compare  the  disk  initial 
instabilities  to  those  of  the  annular  region  in  Bergeron  et  al  [3].  The  annular  region  is 
forced  such  that  the  shear  layer  develops  at  a  distance  of  jR  from  the  center  where  R  is 
the  outer  radius  of  the  annulus,  see  figure  4.1. 


fig  4.1  In  the  annulus  the  shear  layer  is  developed  at  . 

Bergeron  et  al  [3]  define  the  Reynolds  number  using  the  cell  height  for  the  length  scale 
instead  of  the  shear  layer  thickness  so  our  parameters  are  adapted  to  their  definition.  The 
cell  aspect  ratio  is  defined  as 

r-f 

where  a  is  the  distance  from  the  center  of  the  disk  to  the  shear  layer  and  h  is  the  height 
of  the  cell.  We  consider  the  two  extreme  cases  from  Bergeron  et  al  [3]  where  the  aspect 
ratios  are  F  =  6.5  and  F  =  10 ,  respectively.  The  Reynolds  number  defined  in  Bergeron 
et  al  [3]  is 

where  a  is  the  ratio  a  =  for  an  annulus  with  outer  radius  R„  and  inner  radius  R, . 

Regarding  the  disk  geometry  or  =  1  is  the  closest  approximation.  The  normalization  of 
the  disk  sets  R  =  l  and  implies  the  shear  layer  develops  at  a  =  0.75  .  Considering  the 
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r  =  6.5  case  first,  the  resulting  height  is  h=  0.1 154 .  The  viscosity  is  set,  v  =  2  *10  ^ , 
and  now  the  Reynolds  number  is  only  a  function  of  the  angular  velocity,  AQ . 

Applying  forcing  at  the  j  radius  gives  excellent  correspondence  between  the  disk  and 

annulus  regarding  the  onset  of  instability,  see  figure  4.2. 

The  agreement  between  these  simulations  begins  to  deteriorate  at  modes  above 
and  below  the  critical  Fourier  mode,  n  =  5.  Modes  less  than  the  critical  Fourier  mode 
begin  to  exacerbate  the  differences  between  the  disk  and  annular  geometries.  As  shown 
in  section  4.3,  the  disk  and  annulus  dynamics  differ  considerably  for  lower  Fourier  modes 
especially  n  <  2 .  Variations  in  the  modes  above  the  critical  Fourier  mode  are  a  function 
of  increased  angular  velocity  and  the  shear  layer  modeling.  Bergeron  et  al  [3]  use  a  cubic 
spline  to  model  the  velocity  transition  while  we  use  a  hyperbolic  tangent. 

The  second  case  with  F  =  10  provides  even  better  agreement  between  the  two 
stability  analyses.  For  this  second  case  the  height  is  h  =  0.075 .  The  lower  height 
increases  the  Ekman  pumping  coefficient  and  reduces  the  shear  layer  thickness.  The 
marginal  stability  is  shown  in  figure  4.3.  The  variation  in  the  lower  modes  is  consistent 
with  the  previous  diagram.  The  consistency  of  the  homogeneous  instability  demonstrates 
the  closely  related  nature  of  these  simulations  at  low  Reynolds  numbers 

The  experiments  completed  by  Rabaud  &  Couder  [25]  and  Chomaz  et  al  [6]  agree 
with  the  analysis  of  the  initial  instabilities  as  discussed  in  Bergeron  et  al  [3].  The  Landau 
equations,  Bergeron  et  al  [3],  show  that  the  bifurcation  from  steady  flow  is  a  supercritical 
pitchfork  in  the  mean  rotating  reference  frame.  Friih  &  Read  [15]  confirm  the  initial 
instability  as  a  Hopf  bifurcation  in  the  non-rotating  frame.  For  the  remainder  of  the 
results  we  consider  the  shear  layer  at  the  midpoint  of  the  radius,  a  =  0.5 .  The  height  and 
Ekman  pumping  coefficient  are  also  adjusted  to  maintain  an  aspect  ratio  of  ten,  F  =  10 . 
The  complete  list  of  parameters  is  given  here 
;i  =  L  =  0.05 

v  =  2*10"^ 
a  =0.5 
|^  =  6.4 

F  =  10 
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The  marginal  stability  of  the  Fourier  modes  is  shown  in  figure  4.4. 


Marginal  Stability:  Annulus  v.s.  Disk 
-- ->l^ -- Annulus  -o-Disk 


Fourier  Mode 


fig  4.2  Initial  instability  by  Fourier  mode  T  =  6.5  and  a  =  0.75 . 


fig  4.3  Initial  instability  by  Fourier  mode  F  =  10  and  a  =  0.75 . 
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fig  4.4  Marginal  stability  for  a  =  0.5  and  h  =  0.05 . 


The  Reynolds  numbers  used  for  these  calculations  differs  from  the  ones  used  in  Bergeron 
et  al  [3]  per  the  difference  in  definitions.  Recall  we  use  the  Reynolds  number  dependent 
on  the  shear  layer  thickness,  L ,  as  discussed  in  section  4.1 


While  the  actual  values  of  the  Reynolds  number  cannot  be  compared  between  figures  4.3 
and  4.4,  one  notes  the  reduction  in  the  critical  Fourier  mode  number  when  the  shear  layer 
is  located  closer  to  the  center  of  the  disk,  for  a  given  aspect  ratio.  The  same  reduction  in 
critical  mode  number  also  occurs  when  the  aspect  ratio  is  reduced,  Niino  &  Misawa  [23] 
and  Bergeron  et  al  [3].  The  eigenvalue  analysis  described  in  section  3.3  demonstrates  the 
systematic  development  of  the  pseudospectrum  from  right  to  left  as  additional  Chebyshev 
modes  are  retained  in  the  analysis.  Comparison  of  three  resolutions  separates  the 
spurious  eigenvalues  from  those  of  physical  spectrum.  The  resulting  spectrum  for  the 
initial  bifurcation  at  n  =  6  is  shown  in  figure  4.5.  The  spectrum  is  adjusted  to  the  mean 
rotating  frame,  and  the  bifurcation  appears  as  a  supercritical  pitchfork.  The  critical  value 
of  the  Reynolds  number  for  the  onset  of  instability  is  Re  =  37.8 . 

At  Re  =  40  just  beyond  the  initial  bifurcation,  an  n  =  6  braid  begins  to  form.  A 
visualization  of  the  flow  is  provided  in  figure  4.7.  As  the  Reynolds  number  is  increased, 
enstrophy  is  added  to  the  system  increasing  the  strength  of  the  vortices  in  the  braid,  see 
figure  4.17.  The  method  of  increasing  the  Reynolds  number  is  discussed  in  section  4.3. 
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Eigenvalues:  Fourier  mode  6 
384  Chebyshev  modes 


fig  4.5  The  eigenvalues  for  Fourier  mode  6  with  the  pseudospectra  removed. 

The  Ekman  pumping  coefficient  is  a  function  of  height  and  viscosity.  When  a 
shear  layer  remains  constant  and  the  height  is  reduced  to  h  =  0.02 ,  the  marginal  stability 
curve  shifts  due  to  the  increase  in  aspect  ratio  see  figure  4.6.  These  results  are  consistent 
with  the  experimental  work  of  Rabaud  &  Couder  [25]  and  the  analysis  of  Bergeron  et  al 
[3].  For  the  increased  aspect  ratio,  the  critical  mode  moves  from  n  =  6  to  n  =  8  while  the 
critical  Reynolds  number  increases  by  more  than  a  factor  of  three.  Numerical  simulation 
of  the  thin  2D  flow  allows  for  investigation  of  cell  heights  unrealizable  via  an 
experimental  apparatus.  These  investigations  may  lead  to  further  analysis  of  the 
secondary  transitions  for  n  >  8  flows.  Chomaz  et  al  [6]  and  Friih  &  Read  [15]  report 
transitions  at  high  mode  number  which  they  are  unable  to  analyze. 
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Fourier  Mode 

fig  4.6  The  marginal  stability  of  the  rotating  flow  for  h  =  .05  which  gives  the 
Ekman  coefficient  of  ekf  =  6.4  and  h  =  .02  which  gives  ekf  =  25.4 . 


pit/pitl.plf 


fig  4. 7  Flow  at  Re  =  40  just  afier  the  initial  Hopf  bifurcation. 
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4.3  Asymmetric  tripole 

Transitions  of  the  flow  beyond  the  initial  bifurcation  culminate  in  the  asymmetric 
tripole.  The  similarity  between  the  annulus  and  disk  simulations  breaks  down  as  the  flow 
is  driven  to  higher  Reynolds  numbers  and  through  the  secondary  bifurcations.  The  high 
Reynolds  number  flow  leads  to  stability  with  n<2  and  forms  an  asymmetric  tripole. 

The  forms  of  the  tripole  are  only  accessible  within  the  annular  geometry  by  implementing 
artificial  boundary  conditions  that  de-emphasize  the  wall.  Section  4.4  provides  additional 
discussion  on  the  secondary  bifurcations  occurring  prior  to  the  point  vortex. 

The  flow  in  the  disk  geometry  transitions  to  a  single,  coherent  vortex,  “point 
vortex”,  after  the  n  =  3  braid  becomes  unstable.  Experiments  using  annular  regions  have 
unsuccessfully  attempted  to  develop  the  point  vortex  flow.  A  symmetric  balance  exists  in 
the  annular  region,  Niino  &  Misawa  [23],  for  the  n  =  2  solution.  This  balance  also  exists 
for  the  disk  geometry,  Konijnenberg  et  al  [20],  but  as  the  Reynolds  number  is  increased, 
the  symmetric  n  =  2  solution  gives  way  to  the  asymmetric  tripole,  n  =  2  ®  1 .  One  of  the 
outer  vortices  begins  to  collapse  forming  the  point  vortex,  while  the  other  expands  and 
shears.  The  result  is  a  steady  point  vortex  paired  with  a  steady  or  periodic  sheet  like 
vortex  on  the  antipodal  side.  As  the  Reynolds  number  is  increased  the  point  vortex 
continues  to  increase  in  strength.  Figures  4.8  and  4.9  show  a  typical  point  vortex  solution 
at  Re  =  250 .  Note  the  3*^^  Fourier  mode  has  more  energy  than  the  4^**. 

The  parameter  space  for  the  asymmetric  tripole  is  filled  with  complex  dynamics 
from  the  steady  n  =  2  flow  to  the  periodic  n  =  2  ®  1 .  The  rate  of  Reynolds  number 
transition  can  move  the  flow  into  any  one  of  several  attractors.  For  the  numerical 
simulation,  the  Reynolds  number  is  typically  incremented  every  50  iterations  by  a  factor 

of  6.25  *  10”^ .  This  allows  the  flow  50  iterations  to  relax  between  adjustments  while  still 
providing  a  steady  rate  of  change.  A  different  form  of  the  asymmetric  tripole  is  seen  in 
figure  4.10.  This  state  is  found  after  a  period  doubling  bifurcation  from  the  n  =  4  ®  2 
steady  state  to  an  n  =  2  ®  1  steady  state,  figure  4. 1 1.  Here  the  4*  Fourier  mode  has  more 
energy  than  the  3’^'*,  and  one  can  see  the  impact  on  the  form  of  the  tripole.  Although 
beyond  the  scope  of  this  thesis,  our  spectral  simulation  provides  a  tool  for  the  complete 
exploration  of  the  parameter  space  associated  with  the  asymmetric  tripole. 
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fig  4.9  Energy  and  enstrophy  for  the  2  ®  1  point  vortex. 
The  numbers  correspond  to  Fourier  modes. 
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4.4  Secondary  bifurcation 

Transitions  occurring  after  the  initial  bifurcation  offer  a  breadth  of  rich  dynamics. 
The  equilibrium  solution  found  beyond  the  initial  bifurcation  is  a  braid  of  n  vortices.  As 
the  Reynolds  number  is  increased  beyond  the  critical  value,  stability  of  the  initial  vortex 
braid  is  lost,  and  a  transition  to  a  subsequent  vortex  braid  occurs.  With  a  smooth,  slow 
increase  of  the  Reynolds  number  the  new  stable  braid  consists  of  n  - 1  vortices.  Rapid 
transitions  can  bypass  the  n  - 1  attractor  and  lead  to  other  bifurcations  such  as  period 
doubling.  As  the  Reynolds  number  increases,  the  n  - 1  transition  is  often  preceded  by  a 
transition  to  quasi-periodic  states,  and  Chomaz  et  al  [6]  investigate  the  n  =  4  n  =  3 
transition  through  two  symmetry  breaking  bifurcations  and  a  temporal  bifurcation.  Frilh 
&  Read  [15]  perform  a  careful  experimental  study  over  a  wide  range  of  Reynolds 
numbers.  Their  detailed  search  of  the  parameter  space  provides  characterization  of  the 
initial  supercritical  pitchfork  bifurcation,  period  doubling  bifurcations,  irregular  flow 
regions,  and  modal  transitions.  Vortices  with  unusually  high  levels  of  white  noise  are 
reported  for  specific  profiles.  Rand  [26]  analyzes  the  transitions  from  periodic  flows  to 
modulated  (quasi-periodic)  flows.  For  our  purposes  a  periodic  flow  is  defined  as  a  flow 
which  is  periodic  in  time  with  respect  to  the  mean  rotating  coordinate  system.  A  quasi- 
periodic  flow  exhibits  multiple,  spatial  or  temporal  frequencies. 

4.4.1  Period  doubling 

Once  the  flow  is  stable  with  an  even  number  of  vortices,  the  period  doubling 
bifurcation  results  from  the  first  spatial  symmetry  breaking,  Rand  [26],  and  exhibits  a 
direct  exchange  of  energy  between  the  mode  n  and  the  mode  j .  When  the  rate  of 

change  of  Reynolds  number  is  fast  enough,  period  doubling  bifurcations  can  be  found. 
Friih  &  Read  [15]  experimentally  analyze  period  doubling  bifurcations  and  report  smooth 
transitions  during  these  bifurcations.  Our  simulations  confirm  the  smooth  transitions  for 
n  =  6  -4  n  =  3  period  doubling  and  reveal  the  sensitive  nature  of  the  dynamics 

Flow  initially  forced  to  Re  =  75  stabilizes  as  an  n  =  5  braid  and  flow  initially 
forced  to  Re  =  60  stabilizes  as  an  n  =  6  braid.  If  the  fluid  velocity  is  ramped  from 
Re  =  60  to  Re  =  75  at  a  rate  of  1.25/unit  time ,  the  flow  experiences  a  period  doubling 
bifurcation.  This  transition  is  shown  in  figures  4.13  and  4.14.  However,  if  the  flow  is 
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ramped  from  Re  =  60  to  Re  =  67.5  and  allowed  to  stabilize  and  then  ramped  from 
Re  =  67.5  to  Re  =  75 ,  the  period  doubling  does  not  occur  and  the  final  state  is  the  n  =  5 
braid. 

In  order  to  analyze  the  period  doubling  bifurcations,  periodicity  is  enforced  on  the 
flow.  As  an  example  the  flow  at  Re  =  75  is  naturally  n  =  5  braid,  figures  4.15  and  4.16, 
but  can  be  forced  to  maintain  an  n  =  6  braid  by  allowing  only  the  n  =  6k  Fourier  modes 
to  exist,  see  figure  4.17.  The  n  =  6  flow  is  unstable  when  additional  Fourier  modes  are 
considered  in  the  eigenvalue  analysis.  To  introduce  the  spectral  analysis,  we  consider 
stabilized  flow  shortly  after  the  initial  bifurcation.  Re  =  40 ,  figure  4.7.  Here  the 
periodicity  is  not  enforced  since  the  flow  naturally  consists  of  only  n  =  6k  Fourier 
modes.  The  eigenvalue  analysis  for  this  flow  results  in  figure  4.12.  The  flow  is  a  stable 
n  =  6  with  the  leading  order  eigenvalues  having  Re[X\  =  -5.96 .  To  find  bifurcations  to 
quasi-periodic  or  period  doubling  solutions,  we  consider  the  flow  in  an  unstable  region 
with  forced  periodicity. 

Accurate  determination  of  the  instability  requires  inclusion  of  all  the  Fourier 
modes.  This  increases  the  size  of  the  eigenvalue  calculation  from  □  20,000,000 
elements  to  ~  200,000,000  elements.  The  computational  horsepower  required  to  solve 
these  problems  in  not  currently  available.  Therefore,  this  analysis  is  reserved  for  future 
investigation  as  some  additional  shortcuts  may  be  found  to  reduce  the  size  of  the 
eigenvalue  problem  and  computational  capability  is  always  increasing. 
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fig  4.12  Spectrum  for  the  stable  n  =  6  braid  at  Re  =  40 . 
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fig  4.13  Flow  before  and  after  a  period  doubling  bifurcation,  n  =  6  —>  n  =  3 .  The  Reynolds  number  is 
ramped  up  from  60  to  75,  and  all  the  Fourier  modes  are  allowed  to  interact. 


Data  file;  Plot  file: 

17-Nov-O  -  - 

alt/omega.322  plt/pH1.plf 


fig  4.15  Flow  at  Re  =  75  without  forced  periodicity. 
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fig  4.17  Flow  at  Re  =  75  with  periodicity  n  =  6k  enforced. 
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4.4.2  Oscillatory  states 

Transitions  to  periodic  states  has  been  observed  in  various  experiments,  Chomaz 
et  al  [6],  Yao  &  Zabusky  [35],  Friih  &  Read  [15],  and  Konijnenberg  et  al  [20]. 
Konijnenberg  et  al  [20]  and  Yao  &  Zabusky  [35]  report  a  periodic  n  =  3  solution  while 
Chomaz  et  al  [6]  find  the  n  =  4  ®  2  and  «  =  4  ®  2  ®  3  modulated  waves.  In  the  wide 
parameter  search  by  Friih  &  Read  [15],  both  the  periodic  n  =  3  and  n  =  4  ®  2  states  are 
reported.  In  addition,  irregular  flows  are  revealed  for  specific  combinations  of  Rossby 
number  and  Ekman  number.  Near  the  irregular  flow  regimes,  Friih  &  Read  [15]  observe 
“noisy”  vortices  where  "...  the  ‘noisy’  flows  not  only  exhibited  an  increased  level  of 
relatively  uniform  fluctuations,  but  also  frequent  spikes  in  the  signal  of  the  same  order  of 
magnitude  as  the  vortex  signal  itself.”  They  investigated  the  noisy  vortices  by  using 
singular  systems  analysis  (SSA)  to  reconstruct  the  flow  phase  space.  After  performing 
the  SSA  analysis,  they  conclude,  “that  the  noisy  vortices  do  not  contain  any  more 
information  or  further  deterministic  dynamics  (as  could  have  been  indicated  by  the 
enhanced  eigenvalues  of  the  second  pair  of  eigenvectors),  but  that  the  noisy  vortex  state 
contains  the  same  deterministic  dynamics  as  the  periodic  vortices  only  with  added  white 
noise.”  It  is  our  conjecture  that  Friih  &  Read  [15]  were  near  the  n  =  4  ®  2  ®  3  periodic 
state,  but  were  not  close  enough  to  isolate  the  mode  3  contribution. 

The  n  =  4®2®3  is  seen  by  Chomaz  et  al  [6]  as  the  second  symmetry  breaking 
bifurcation  when  the  flow  transitions  from  n  =  4  — >  n  =  3 .  We  have  found  the  state  to 
exist  as  a  periodic  solution  for  a  small  range  of  Reynolds  numbers.  The  time  evolution  of 
the  n  =  4  ®  2  and  n  =  4  ®  2  ®  3  periodic  states  is  shown  in  figures  4.19  and  4.21, 
respectively.  The  associated  energy  and  enstrophy  evolutions  are  included  in  figures  4.20 
and  4.21.  All  transitions  are  assumed  to  be  complete  since  the  vortex  structures  have  not 
changed  for  over  250  cycles  of  their  respective  turn  around  times.  In  addition  to  the 
reports  by  Chomaz  et  al  [6]  and  Friih  &  Read  [15]  we  find  two  distinct  n  =  4  ®  2  periodic 
states.  One  occurs  before  the  first  symmetry  breaking  bifurcation  where  the  mode  2 
appears,  while  the  other  occurs  after  the  first  symmetry  breaking  bifurcation  and  prior  to 
the  second.  Chomaz  et  al  [6]  also  reports  a  transient  “chaotic”  state  analogous  to  the 
irregular  flows  seen  by  Friih  &  Read  [15].  We  have  observed  the  irregular  flow  at  a 
Reynolds  number  slightly  higher  than  the  n  =  4  ®  2  ®  3  state,  figures  4.25  and  4.26. 
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These  irregular  flows  may  be  either  quasi-periodic  in  the  rotating  frame  or  chaotic.  Time 
series  analysis  of  the  data  is  required  to  determine  and  dimensionality  of  the  dynamics 
and  is  beyond  the  scope  of  this  thesis. 


Thus,  we  provide  numerical  simulation  for  the  full  series  of  bifurcations  between 
the  n  =  4  steady  state  and  the  n  =  3  steady  state.  Figure  4.18  is  a  cartoon  of  the 
transitions  between  the  steady,  periodic,  and  irregular  states. 


fig  4.18  These  are  the  states  surrounding  the  4  ©  3  ©  2  periodic  solution. 

The  periodicity  is  time  dependent  with  s  steady  and  p  =>  periodic . 
Unverified  results  have  been  reported  in  other  experiments  and  simulations. 


The  bold  lines  indicate  levels  at  which  hysteresis  is  observed  e.g.  once  a  steady  or 
periodic  n-3  flow  is  established  the  symmetry  breaking  bifurcations  between  the  n  =  4 
and  M  =  3  braids  cannot  be  reversed  by  lowering  the  Reynolds  number.  The  transitions 
consist  of  spatial  and/or  temporal  symmetry  breaking. 

With  the  highly  accurate  spectral  simulation  additional  investigations  can  be  made 
into  the  bifurcations  between  the  n  =  3  and  n  =  2  e  1  states.  In  order  to  isolate  the 
bifurcations  and  associated  states,  the  flow  must  be  allowed  to  relax.  All  intermediate 
periodic  states  are  bypassed  if  the  flow  does  not  relax  for  at  least  20  cycles  of  the  turn 
around  time  and  some  bifurcations  require  more  than  200  cycles. 


alt/omega.351  pll/pmz.ptf 


fig  4.19  4  ®  2  periodic  flow. 
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fig  4.20  Energy  and  enstrophy  for  the  4*2  periodic  state. 
The  numbers  correspond  to  Fourier  modes. 
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fig  4.23  Asymmetric  fiow  resulting  from  the  4  ©  3  ®  2  periodic  state. 
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fig  4.24  Energy  and  enstrophy  for  the  transition  from  4  ®  2  steady  to  the  4  ©  3  ©  2 
periodic.  The  numbers  correspond  to  the  Fourier  modes. 


t/omega.970 


fig  4.25  Irregular  fiow  after  the  4  ®  3  ®  2  periodic  state. 
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APPENDIX  A 


Vector  identities 

A.l  ax{b'x.c^  =  {a-c^b -{a-b^c 

A.2  •  6  ^  =  a  •  V6  +  5  •  V5  +  a  X  X  5)  +  6  X  (V  X  a) 

A.3  V^a  =  V(V-a)-Vx(Vx5) 

A.4  Vx(a5^  =  Vax6+a(Vx5) 

A.5  V •{ab^  =  a^ -b +  b -Va 

A.6  Vx^Sx5)  =  5-V5-6(V-a)-a-V6 

A.7  V-(5x5)  =  6-(Vxa)-a-(Vx5) 
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APPENDIX  B 


Conditioning  of  the  problem 

In  general  the  linear  Helmholtz  equation  has  the  form 
B.l  A^”=| 


where  a  is  used  to  find  the  vorticity  through  the  smoothing  and  postconditioning 
B.2  ®  =  for  7=1,2 

To  estimate  the  condition  of  the  problem  consider  a  perturbed  linear  system  B.l  and  a 
perturbed  conversion  B.2  giving 

B.3  {A  +  SA)[3°  +S3'‘)  =  g  +  Sg 

B.4  &  +  Sa  =  (r^  +  SR^ )  (P  +  <^P)  (d)°  +  ) 

The  conditioning  for  the  problem  is  measured  through  the  magnitude  of  the  relative  error 
in  vorticity.  This  error  compares  the  norm  of  the  vorticity  error  to  the  norm  of  the 
vorticity. 

For  the  first  order  approximation  all  quadratic  error  terms  are  eliminated.  To 
obtain  the  condition  number  for  the  Helmholtz  operator,  start  by  solving  B.3  for  S3°  and 
,  respectively 

B.5  33°  =  A~^g  +  A~^Sg  -3°  -  A~^SA3°  =  A~^Sg  -  A~^SA3° 

B.6  3°  =A"'^{g  +  Sg)-S3° 


Assuming  A  ^  ||(5^A||<1,  recursive  application  of  B.6  into  B.5  while  eliminating  higher 
order  terms  gives 


B.7 


33°  =  A  ^3g  + 


.k=l 


A  {g-^Sg) 


Now  expand  B.4  while  subtracting  the  representation  for  vorticity,  3  =  R^P^®",  giving 
33  =  (r^<5P^  +  (5R^P^ )»  +  (r^P^  +  R^<?P^  +  3R^Y^)33° 

Substituting  B.7  into  this  equation  gives 
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da  =  (r^'<5P^  +  )  a  +  (r^P^  +  +  SR^?^ )  A"'^^  + 


B.8 


(r^P^  +  R^^P^  +  ^R^P^ )  X(-l)*  (a‘*^a)* 

Vk^l 


A"’  (l  +  ^g) 


The  summation  is  transformed  in  two  steps  enabling  an  efficient  norm.  Note  all  elements 

in  the  sum  have  a  A~^  product  on  the  left  and  a  ^A  product  on  the  right.  Factoring  out 
these  common  terms  gives 

B.9  |;(-1)‘  (a-Va)‘  =  A-'  [  (5AA-'  f 

k=\  La=1 

The  sum  on  the  right  hand  side  is  now  in  geometric  form.  When  a  norm  is  applied  to  the 
geometric  sum  an  inequality  follows 


^A 


|;(-i)*(<jAA-)‘- 

*=1 


1 


\-5AA- 


1 


i-IMI 


When  B.9  replaces  the  sum  in  B.8  and  a  norm  is  taken  from  both  sides,  the  result  is  an 
inequality  of  the  form 


B.IO 


IMN 

R^  l  JP^ 

+  ||^R^  P^  ( 

r  +1 

R^'P^A"’  1 

Ir-' 

||^p2 

A” 

‘  ||5g||+||5Rf  ||P" 

A-' 

INF 

(||r-'p^a"‘ 

+ 

R^ 

A"‘||  + 

P"  A-‘  ) 

— o 

G) 


+  A- 


Il'Jfll) 


where  S  = 


-4 — Til .  At  this  point  we  introduce  a  notation  to  simplify  the  above 

estimate.  The  Chi  function  is  the  ratio  of  the  norm  of  the  error  to  the  norm  of  the 
operator  e.g.  for  an  operator  ^  the  Chi  function  is  defined  as 

y(P\  =  JIM 

1^1 


Applying  the  Chi  function  while  dividing  by 


CD 


transforms  B.IO  into 
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II5< 


|||p1|[;ir(p^)H-z(R0] 


m 


B.ll 


I!  II 

a-IIIr-'IIIIp 
ISIIII^AI 

A 


+ 


ll^g|l  IIp7T>2  a  -1 


R^P^A"M  + 


1[z(p')+z(rO/ 


r^p^a"H11  1  + 


A-‘ 

INI) 

«  ) 

R^ 


ii[^(p^)+z(RO]( 


1+ 


A  -ill  IN 

m. 


Further  simplification  comes  through  the  norm  inequality  based  on  B.l 

_!_<M 

iKrilfii 

and  an  additional  definition,  the  condition  number  of  A , 


Ca=I|a 


-1 


Rewriting  B.l  1  with  these  considerations  gives 


IMI 


R^ 


i|[z(P')+z(R0‘- 


+  IIAII  R-'P^A”* 


\\zig)- 


1r1I11p1|z  (g)[z  (p" ) + z  (r' )] + 

1A|R|Z(A)1R^P"A-'||[1  +  C^;ir(g)]  + 

||Hz(A)[z(P^)+z(R')][i+CAz(g)] 


CaP||||p 

For  the  final  estimate  remove  the  high  order  terms  and  neglect  x  )  ^nd  %  )  • 

Since  the  matrices  R^  and  P^  have  error  controlled  by  the  truncation  of  the  Chebyshev 
expansion,  the  errors  are  exponentially  small  when  compared  to  those  of  the  SBDF 
scheme  found  in  A .  These  simplifications  lead  to 

M<||A||||r^P^A-|;,(|)+||A||12|;,(A)||r^P^A-' 


B.12 


||A||  R'P^A-‘  [;tr(g)+ll|^(A)] 


The  ratio  of  norms  between  the  calculated  value,  a>  ,  and  the  vorticity  forms  the 
inequality 


and  this  inequality  applied  to  B.12  giving  the  final  form  of  the  condition  number  bound 
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M<ll|l||R>P^A-‘||[;f(|)  +  lSl/(A)]<M||R'P^A-‘||[z(l)  +  2z 

This  bound  is  nearly  the  product  of  the  condition  number  of  A  with  the  relative  error  in 
the  problem. 
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APPENDIX  C 


Chebyshev  recurrence  relations 

The  Chebyshev  polynomials  provide  an  orthogonal  basis  for  fast  computation  in 
numerical  simulation.  In  addition  to  possessing  a  fast  transform,  Fourier  Cosine 
Transform  (FCT),  these  polynomials  relate  through  several  recurrence  relations.  The 
Chebyshev  polynomials  are  defined  by  the  relation 
C.1  7;(r)  =  2rr,.,(r)-r..,(r) 

where  Tq=  \  and  Ty{x)  =  x .  An  alternative  identity  for  the  polynomials  relates  to  their 
cosine  transform 

C.2  (cos  =  cos  mO 

Considering  the  Chebyshev  polynomials  within  the  range  -1  <  r  <  1  provides  additional 
relations.  At  the  boundaries  the  polynomials  are  unitary 

r,(-l)  =  (-!)" 

and  at  the  middle  of  the  interval,  r  =  0 ,  the  behavior  is 
r  0  m  odd 
m  =  2i 

The  zeros  of  the  Chebyshev  polynomials  are  distinct,  and  (r)  has  exactly  n  zeros  in 
-1  <  r  <  1 .  The  form  for  the  zeros  is 
^  =  0  for  rjt  =  cos(^;r) 
where  k  =  0,l,...,n  - 1 . 

Multiplying  two  Chebyshev  polynomials  produces  a  relation  through  C.2. 
Consider  the  product  of  polynomials 

T„  (cosd)T^  [cosO)  =  cosnOcosmO 
and  using 

cosnOcosmO  =  j[cos(n^  +  mO)  +  cos(^n0  - 
gives  the  relation 
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The  product  relation  is  used  to  derive  the  Chebyshev  convolution  in  appendix  E. 

Differentiation  of  a  Chebyshev  polynomials  also  produces  a  recurrence.  The 
derivative  of  C.l  with  respect  to  r  gives 

recursive  substitutions  of  this  form  while  shifting  the  index  give 

r  =m(2T  .  +%) 
m  y  m-i  m-2  J 

where  T(=Tq  and  Tq  =  0 .  The  differentiation  recurrence  produces  the  Chebyshev 
differentiation  matrix. 
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APPENDIX  D 


Chebyshev  matrices 

Multiplication  and  differentiation  by  the  radius,  r ,  appear  throughout  the  Navier- 
Stokes  equations.  When  the  Fourier-Chebyshev  expansions  are  truncated  the  problem 
can  be  solved  using  matrix-vector  representation.  Flow  parameters  are  represented  by 
their  Fourier-Chebyshev  coefficients.  Since  the  problem  imcouples  through  the  Fourier 
modes,  Chebyshev  coefficients  for  a  specific  Fourier  mode  are  rendered  in  a  single  vector 

=  C?,l  Cj,2  ■■■  Ci,m] 


where  is  the  vector  of  Chebyshev  coefficients  for  Fourier  mode  n .  The  required 

linear  operators  in  the  Chebyshev  basis  are  represented  as  matrices.  Multiplication  by  r 
is  nothing  more  than  convolution  with  the  coefficient  vector 

r  =  [0  1  0  •••  of 

Chebyshev  convolution  matrices  are  discussed  in  appendix  E.  The  convolution  matrix 
for  r  is 


R  = 


0  i 
1  0 

0  i 
0  0 
0  0 


0 

1 

2 


0  0 
0  0 

0  ■•.  0 

•.  •.  i 

•  2 

0  4-0 


Differentiation  by  r  derives  from  the  recurrence  relation 


where  T{=Tq  and  Tq  =  0 .  Expanding  the  recurrence  and  truncating  at  mode  M  gives 


0 

1 

0 

3 

0 

5  ••• 

M 

0 

0 

4 

0 

8 

0 

0 

0 

0 

0 

6 

0 

10 

2M 

0 

0 

0 

0 

8 

0 

0 

0 

0 

0 

0 

0 

10 

2M 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2M 

0 

0 

0 

0 

0 

0  0 

0 
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Both  R  and  D  apply  to  a  single  Chebyshev  coefficient  vector.  When  the  full  vector 
form  is  required  as  with  the  secondary  stability  analysis,  composition  of  matrices  is 
required.  The  full  vector  form  of  the  flow  parameters  follows 

fo'W 

0 

;’('•) 

?!('■) 

where  ^  ij',0^  is  either  vorticity,  a ,  or  stream  function,  ^ ,  and  are  the 

Chebyshev  coefficient  vectors  for  the  Fourier  mode,  cosine  and  sine,  respectively. 
Since  multiplication  and  differentiation  by  r  impact  only  the  Chebyshev  coefficients,  the 
matrix  providing  multiplication  by  r  for  the  long  vector  problem  is 

■r  0  •••  0 

0  R  i 
R  =  I^+i®R=  .  . 

0  •••  0  R 

Each  R  matrix  has  size  (M  + 1)  x  {M  + 1) ,  and  there  are  iV  +  2  R  matrices  along  the 

diagonal.  Similarly  differentiation  by  r  produces  the  matrix 

"D  0  •••  0 
0  D 

»  =  V,®“=  ;  0 

0  0  D 

Each  D  matrix  has  size  (M  + 1)  x  {M  + 1) ,  and  there  are  iV  +  2  D  matrices  along  the 
diagonal.  In  addition  to  these  basic  operations  within  the  Fourier-Chebyshev  expansion, 

A  A 

R  s  r  and  D  s  5,. ,  Fourier  modes  require  differentiation  by  ^ .  In  the  long  vector 


environment  a  new  matrix  is  defined  for  this  purpose,  Q  =  0^ .  The  Q  matrix  is  formed 

by  considering  the  differentiation  of  sine  and  cosine  functions.  Here  the  impact  is  strictly 

100 


upon  the  Fourier  modes  leaving  the  Chebyshev  expansions  unchanged.  For 
differentiation  by  9 


where  n-\  equates  to  the  row  of  I^r 
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APPENDIX  E 


Chebyshev  convolution  matrix 

The  Chebyshev  convolution  matrix  is  used  with  the  stability  analysis.  The 
convolution  matrix  is  built  using  the  Chebyshev  multiplication  recurrence  relation 


Additional  notation  is  now  introduced  to  provide  a  matrix-vector  representation  for  the 
Chebyshev  polynomials.  After  the  expansions  are  truncated,  every  polynomial,  Tj ,  can  be 

written  as  a  matrix,  ,  which  is  defined  as: 

<r-  7*  + 1  row 


Here  the  convolution  of  two  Chebyshev  vectors  is  done  only  once  so  a  computationally 
expensive  matrix-vector  multiplication  is  used.  In  general  the  convolution  of  two 
Chebyshev  vectors  g  and  h  is  accomplished  in  the  following  manner 

where  each  Tjh  is  a  matrix  vector  multiplication,  and  gj  is  the  coefficient  for  the  f' 
Chebyshev  polynomial  matrix,  Ty . 

An  example  for  our  problem:  let 

=  [^«,0  ^n,\  ^n,2 

be  the  coefficient  vector  then  the  Chebyshev  convolution  matrix  for  ^  is  given  by  the 
following  matrix 


T,= 


0  1 
1  0 


0  i 


0^0 


1 
2 

0  4-0 


0 

0  4-0 


0  1 


0  i 


0 


0 


Note  due  to  the  non  negative  entries  along  the  positive  slope  diagonal  in  each  matrix, 

_  M 

the  representation  d^co„  Q  =  is  only  accurate  for  the  0  ^  ^  Chebyshev  modes. 

y=o 

Convolution  matrices  created  using  0  ->  M  modes  and  larger  than  (^  + 1)  x  (-y  + 1)  are 

missing  meaningful  data  contributed  by  the  M  +  \  and  higher  modes.  Thus  an  M  mode 
accurate  Chebyshev  convolution  matrix  requires  2M  modes  of  data  within  the 

coefficient  vector. 
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APPENDIX  F 

Fourier  convolution  matrix 

As  with  the  Chebyshev  basis,  the  Fourier  convolution  of  two  functions  generates 
linear  operator  by  representing  a  vector  of  Fourier  coefficients  as  an  equivalent  matrix. 
The  convolution  of  two  vectors  within  the  Fourier  trigonometric  basis  is  more  easily  seen 
by  converting  the  convolution  represented  in  the  exponential  basis.  Therefore,  the 
Fourier  convolution  is  initially  derived  for  an  exponential  basis. 

The  Fourier  convolution  matrix  resulting  from  an  exponential  basis  rests  upon  the 
product  rule  for  exponential  expressions.  Consider  two  expansions  represented  by 
exponential  bases 


N 


The  product  of  these  two  functions,  p  =  cy,  can  be  represented  within  the  same 
exponential  basis 


Every  coefficient,  ,  is  a  sum  of  coefficient  products  e.g. 


Po = c.KriL + +  •  •  • + co/o  +  • " + + c^r-ii 

Using  matrix- vector  representation  all  of  the  functions  are  written  as  vectors  containing 
their  respective  coefficients,  p ,  c ,  and  f .  The  functions  c  and  y  in  vector  form 


- 

-|T 

c  = 

C  N 

^  2 

2 

^1 

11 

1  1 

1 _ 1 

lead  to  the  following  convolution  matrix  and  product 
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’  ^0 

C-1 

C_2 

... 

C  AT 

2 

0 

0 

... 

0 

Cl 

Co 

^-f+1 

C  N 

2 

. 

0 

Cf-1 

... 

... 

Co 

... 

... 

... 

2 

0 

r-i 

Cjn 

2 

Cm  1 

2  ^ 

Cl 

^0 

C-l 

C-f+l 

C  N 

2 

Vo 

0 

2 

Co 

C-f+l 

n 

0 

... 

0 

Cn 

2 

Cn  1 

2  ^ 

... 

Cl 

Co 

C-l 

0 

... 

... 

0 

Cn 

2 

... 

C2 

Cl 

Co 

.  . 

The  result  of  this  matrix  vector  multiplication  is  the  vector  p .  Each  coefficient  is  a 
vector  dot  product  given  by 

IN+n 

Pn=Y.^*n-krk-iL 

F.2 

2N-n 

=  E  ^f-krk-f^„  0>n>f 

k=0 

It  is  important  to  note  that  the  only  coefficient  using  all  of  the  available  information 
from  both  c  and  f  is  Pq.  The  other  coefficients  in  p  are  less  accurate  due  to  the 
truncated  expansions  of  c  and  f .  Thus  to  have  an  accurate  Fourier  convolution  matrix 
for  modes  -  y  <  «  <  y  requires  a  product  of  series  having  size  -N  <n<N .  As  with  the 

Chebyshev  convolution,  the  Fourier  convolution  matrix  requires  a  vector  twice  the  size  of 
the  final  matrix.  The  actual  convolution  matrix  representing  c  is 


105 


^0 

C-1 

C-2 

... 

2 

C_lL-.\ 

2  ^ 

... 

C-V 

'  ^-4  ' 

P  N 

Cl 

Cq 

C-1 

... 

C-f+l 

C  N  C  N  \ 

2  2  ‘ 

... 

C-Af+1 

2 

2 

A’-f+i 

Cf-1 

... 

... 

... 

... 

C  N 
“*  2 

C-f-1 

r-i 

P-x 

Cji 

2 

Cf-1 

Co 

C-f+l 

C_1L 

2 

ro 

= 

Pq 

Cf+l 

Cn 

2 

... 

^  ... 

... 

^0 

... 

C-f+l 

Yx 

Px 

CjV-1 

.  .  . 

Cf+1 

Cn 

2 

Cn  1 

2  ^ 

Co 

C-l 

YiL-x 

A’f-i 

Cv 

... 

... 

Cal+1 

Cn 

2 

C2 

Cl 

Co 

.  . 

pit 

where  the  result  is 

F-3  P„  =  for  -f  <  n  < f . 

k=0 

With  the  convolution  matrix  defined  for  an  exponential  basis,  we  need  convert 
from  the  exponential  basis  to  the  trigonometric  basis.  Using  the  standard  form  for  the 
Fourier  trigonometric  expansion 

N  N 

F.4  ^  =  ^0  +  ^  ^os  nd  +  bf,  sin  nO 

n=-N  n=l 

provides  three  cases  to  investigate 

1)  «  =  0  —  resulting  in  =  Cq 

N  N 

2)  l>n>N  —  then  (oosn^  +  /sin«^) 

n=\  «=i 


3)  -N  <n<-l  (but  make  the  following  change  of  variables:  n  =  -n  so  the  sum 

N  N 

is  1  ^  as  in  case  2)  -  then  ^c_„e"'"'^  =  ^c_„  (cosnO  -isinnO) 


n=l 


n=l 


Combining  the  above  cases  leads  to  the  expansion 

N  _  N 

F.5  Y,  =  ^0  +  Z[(^«  +  c_„)cosn6»  +  i{c„  - c_„)sin«^] 

n=-N  n=\ 

The  results  of  interest  obtained  by  equating  F.4  and  F.5  are  +  ^-n » 

b„  =  i (c„  Rewriting  the  coefficient  equations  while  solving  for  c„  and  c_„  gives 
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exponential  and  trigonometric  Fourier  series  in  combination  with  the  exponential  Fourier 
convolution  matrix  gives  the  following  trigonometric  Fourier  convolution  matrix  and 
product 
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The  trigonometric  Fourier  convolution  matrix  consists  of  four  distinct  blocks  with  each 
block  having  size  ( y  +  l)x»  +  l) .  Each  blocks  is  produced  by  the  cosine  or  sine 


107 


coefficients  of  a  Fourier  vector  containing  + 1  entries.  The  Fourier  convolution  matrix 
applies  to  vectors  written  in  the  full  form.  Consider  the  generic  expansion 

K 

C (r)cosn^  (^)sinn^) 

n=0 

where  ^  is  either  stream  function  or  vorticity  or  »)  and  (r)  and  (r)  are  the 

Chebyshev  coefficient  vectors  for  the  Fourier  cosine  and  sine  modes.  Translating  the 
above  expression  into  a  veetor  of  coefficients  yields  the  full  vector  form 

f(''.»)=[fo('')  -  fiw  0  ?;(>■)  flW  - 

where  the  elements  for  the  Fourier  vector  are  the  Chebyshev  coeffieient  vectors.  To 
produce  a  convolution  matrix  for  the  full  vector  requires  two  steps.  First  the  Fourier 
convolution  matrix  is  constructed  from  the  Chebyshev  coefficient  vectors.  As  a  result 
every  element  of  the  Fourier  convolution  matrix  is  a  coefficient  vector.  In  the  second 
step  each  Chebyshev  coefficient  vector  is  expanded  as  a  convolution  matrix  per 
appendix  E. 

An  example  for  our  problem:  consider  the  expansions  for  the  perturbed  vorticity 
and  stream  function 


)  0  •••  <S|.i(r)] 

l^l(^^)  =  [l^0;l('')  1^1-1  W 

)  ••• 

)  0  y/l^{r)  ••• 

The  convolution  terms  are  entirely  represented  by  considering  *Dg>,  .  This  term  is 

a  convolution  of  the  vectors  RQ?i?o  •  The  application  of  the  Q  matrix  swaps  the 

cosine  and  sine  terms  in  while  applying  the  appropriate  scaling  and  sign  adjustment 

Ql^O=[0  •••  0  -'^1,0  -21^2;0  •" 

A  A 

Since  the  R  and  D  matriees  are  block  diagonal  they  apply  to  the  irmer  Chebyshev 
vectors  without  introducing  any  cross  vector  interaction  within  the  Fourier  space.  The 
final  forms  of  the  RQ^/q  D<y,  vectors  are 

RQ|(;o=[o  R#,'#  2R!<^2'o  -  -2R!i^2'o  - 
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D«,=|^D<yo.i  D6)2;i  •"  D®|.,  0  D«2;i  •" 

Returning  to  the  convolution,  RQ{«?o  *D®, ,  the  Fourier  convolution  matrix  for  RQ^i^o 
be  constructed  in  the  manner  described  above.  This  convolution  matrix  is  written  as 
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and  the  blocks  of  this  matrix  are 
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However,  this  representation  is  misleading  as  each  element  of  F  appears  to  be  an 
(a/  + 1)  X 1  Chebyshev  vector.  In  reality  these  vectors  must  be  replaced  by  the 


appropriate  (M  +  1)x(M  +  1)  Chebyshev  convolution  matrices  e.g.  must  be 

replaced  by  the  Chebyshev  convolution  matrix  R(^f;o  using  the  method  described  in 
appendix  E.  Recall  the  definitions  of  the  Fourier  vector  functions 
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M 

/w=0 

M 

W=:0 

The  Chebyshev  convolution  matrix  resulting  from  the  vector  is  given  by 

IM 

\.,=Zw:«o- 

m=0 

The  Y_c^  matrix  is  cropped  to  (M  + 1)  x  (M  + 1)  as  to  contain  the  full  data  for  a 
(2M  +  l)xl  vector.  has  a  similar  convolution  matrix,  . 


After  the  Fourier  convolution  matrix,  RQ^^o » formed,  multiplication  by  the 
Chebyshev  differentiation  matrix  yields  RQ^^oD.  Since  D  is  block  diagonal  each 
element  of  RQ^gD  has  the  form  Y_c  D  or  Y.,  D ,  a  matrix  product  of  the  Chebyshev 
convolution  matrix,  Y.c  or  Y_,  ,  and  the  Chebyshev  differentiation  matrix  D .  This  is 
the  final  form  of  the  matrix  as  a  product  with  the  unknown  vector,  d)i . 


no 
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